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ABSTRACT 


Linear least squares estimation and regression analyses continue 
to play a major role in orbit determination and related areas. In 
this report we document a library of FORTRAN subroutines that have 
been developed to facilitate analyses of a variety of parameter 
estimation problems. Our purpose is to present an easy to use multi- 
purpose set of algorithms that are reasonably efficient and which use 
a minimal amount of computer storage. Subroutine inputs, outputs, 
usage and listings are given, along with examples of how these routines 
can be used. The following outline indicates the scope of this report; 
Section I, introduction with reference to background material; Section 
II, examples and applications; Section III, a subroutine directory 
summary; Section IV, the subroutine directory user description with 
input, output and usage explained; and Section V, subroutine FORTRAN 
listings. The routines are compact and efficient and are far superior 
to the normal equation data processing algorithms that are often used 
for least squares analyses. 
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I. Introduction 

Techniques related to least squares parameter estimation play a 
prominent role in orbit determination and related analyses. Numerical 
and algorithmic aspects of least squares computation are documented 
in the excellent reference work by Lawson and Hanson, Ref. [1], Their 
algorithms, available from the JPL subroutine library. Ref. [2], are 
very reliable and general. Experience has, however, shown that in 
reasonably well posed problems one can streamline the least squares 
algorithm codes and reduce both storage and computer times. In this 
report, we document a collection of subroutines most of which we have 
written that can be used to solve a variety of parameter estimation 
problems . 


The algorithms for the most part involve triangular and/or 
symmetric matrices and to reduce storage requirements these are stored 
in vector form, e.g., an upper triangular matrix U is written as 


Uii U 

u, 


12 

"l3 

"l4 

^22 

"23 

"24 


"33 

"34 



"44 


etc. 


U(l) 



U(2) 

U(4) 

U(7) 

U(3) 

U(5) 

^Wetc. 


U(6) 

U(9) 



U(10) 


Thus, the element from row i and column j of U, i < j> is stored in 
vector component j(j“l)/2 + i. We hasten, to point out that the engineer, 
with few exceptions, need have no direct contact with the vector sub- 
scripting. By this we mean that the vector subscript related operations 
are internal to the subroutines, vector arrays transmitted from one 
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subroutine to another are compatible, and vector arrays dxsplayed 
using the print subroutine TRIMA.T appear in a triangular matrix format. 

Aside ; The most notable exception is that matrix problems are generally 

formulated using doubly subscripted arrays. Transforming a double 

subscripted symmetric or upper triangular matrix A(-j-) to a vector 

stored form, U{*) is quite simply accomplished in FORTRAN via 
IJ = 0 

DO 1 J = 1,N 
DO 1 I = 1,J 
IJ - IJ+1 
1 ' U(IJ) = A(I,J) 

Similarly, transforming an initial vector D(*) of diagonal positions of 
a vector stored form, U(")> is accomplished using 


JJ = 0 

DO 1 J = 1,N 
JJ = JJ+J 
1 U(JJ) = D(J) 


JJ = N*(N+l)/2 
DO 1 J = N,l,-1 
U(JJ) = D(J) 

1 JJ = JJ-J 


The conversion on the right has the modest advantage that D and U can 
share common storage (i. e. , U can overwrite D) . These conversions 
are too brief to be efficiently used as subroutines. It seems that when 
such conversions are needed one can readily include them as in line code. 
End of Aside 

r 

Although this package of subroutines is designed in the main, for 

✓ 

the analysis of parameter estimation problems one can use it to solve 
problems that involve process noise. With modest amounts of additional 
programming one can even apply our package to filtering problems that 
involve colored noise and mapping. In the latter case, however, reduc- 
tions gained from our use of vector storage are for the most part lost . 
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Mathematical background regardxng Householder orthogonal trans- 
formations for least squares analyses and U-D matrix factorization 
for covariance matrix analyses are discussed in references [1] and [3]. 
Our plan is to illustrate, in Section II, with examples how one can 
use the basic algorithms and matrix manipulation to solve a variety 
of important problems. The subroutines which comprise our estimation 
subroutine package are described in Section III, and detailed input/ 
output descriptions are presented in Section IV. 

Section V contains FORTRAN listings of the subroutines. There are 
several reasons for including such listings. Making these listings 
available to the engineer analyst allows him to assess algorithm 
complexity for himself; and to appreciate the simplicity of the 
routines he tends otherwise to use as a black box. The routines are 
not truly portable, and users can, when necessary make modifications 
so that the subroutine package can operate on systems other than the 
UNIVAC 1108. When estimation problems arise to which our package does 
not directly apply (or which can be made to apply by an awkward conca- 
tenation of the routines) one may be able to modify the codes and widen 
still further the class of problems that can be efficiently solved. 
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II. APPLICATIONS AND EXAMPLES 

Our purpose in this section is to illustrate, with a number of 
examples, some of the problems that can be solved using this ESP. The 
examples, in addition, serve to catalogue certain estimation techniques that 
are quite useful. 

To begin, let us catalogue the subroutines that comprise the ESP: 


1) 

AGTSN 

(A G Turner) 

Agee-Turner rank 1 update 

2) 

A2A1 

(A to A one) 

Matrix A to matrix A1 

3) 

COMBO 

(combo) 

Combine R and A namelists 

4) 

C0V2RI 

(cov to R I) 

Covariance to R inverse 

5) 

C0V2UII 

(cov to U D) 

Covariance to U-D factors 

6) 

C2C 

(C to C) 

Permute the rows and columns of matrix C 

7) 

INF2R 

(inf to R) 

Infoinnation matrix to (triangular) R 

8) 

PERMUT 

(permute) 

Permute the columns of matrix A 

9) 

RINCON 

(rin con) 

R inverse with condition number bound 

10) 

RI2C0V 

(R I to cov) 

R inverse to covariance 

11) 

R2A 

(R to A) 

Triangular R to matrix A 

12) 

R2RA 

(R to R A) 

Transfer a triangular block of R to trian- 
gular RA 

13) 

RUDR 

(rudder) 

SRIF R to U-D factors or vice versa 

14) 

THH 

(T H H) 

Triangular Householder data processing 

15) 

TRIMAT 

(trl mat) 

Triangular matrix print 

16) 

TTHH 

(TTHH) 

Two triangular matrix Householder processing 

17) 

TZERO 

(T zero) 

Zero a horizontal segment of a triangular 
matrix 
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18) 

UDMES 

(U D measurement) 

U-D measurement updating 

19) 

UD2C0V 

(U D to cov) 

U-D factors to covariance 

20) 

UD2S1G 

(U D to sig) 

U-D factors to sigmas 

21) 

UTINV 

(U T inverse) 

Upper triangular matrix inverse 

22) 

UTIROW 


Upper triangular inverse, inverting only 
the upper rows 

23) 

WGS 

(W G-S) 

Weighted Gram-Schmidt triangular reduction 


These routines are described in succeedingly more detail in sections III, 
IV, and V. The examples to follow are chosen to demonstrate how these 
various subroutines can be used to solve orbit determination and other 
parameter estimation problems. It is important to keep in mind that these 
examples are not by any means all inclusive, and that this package of 
subroutines has 'a wide scope of applicability. 

II. 1 Simple Least Squares 

Given data in the form of an overdetermined systems of linear 
equations one may want a) the least squares solution; b) the estimate 
error covariance, assuming that the data has normalized errors; and 
c) the sum of squares of the residuals. The solution to this problem, 
using the ESP can be symbolically depicted as 

• [A z] ►[R z] , e 

Remarks : The array [A z] corresponds to the equations Ax = z-v, veN(0,I); 

the array [R z] corresponds to the triangular data equation Ex = z-v, 
vsN(0,I) and e = |)z-Ax|j 


[R z] 


UTINV 


[R x] 


Remark ; 



z 
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One may be concerned with the integrity of the computed inyerse 
and the estimate. If one uses subroutine RINCON instead of UTINy then 

in addition one obtains an estimate (lower and upper bounds) for the 

condition number R, If this condition number estimate is large the 

computed inverse and estimate are to he regarded with suspicion. By 

large, we mean considerable with the machine accuracy (ylz. on an 18 

15 

decimal digit machine numbers larger than 10 ) . Note that the condi- 

tion number estimate is obtained with negligible additional computation 
and storage. 




Remarks ; C - R R = estimate error covariance. Some computation can 
be avoided in RI2COV if only some (or all) of the standard deviations 
are wanted. 

II. 2 Least Squares With A Priori 

If a priori information is given, it can be included as additional 
equations (in the A array) or used to initialize the R array in subroutine 
THH (see the subroutine argument description given in section IV) . One is 
sometimes interested in seeing how the estimate and/or the formal 
statistics change corresponding to the use of different a priori 
conditions. In this case one should compute [R z] as in case II. 1, and 
then include the a priori [R^ z^] using either subroutine THH, or 
subroutine TTHH when the a priori is diagonal or triangular, e.g. , 




[R 

[R z }) 
o o' 


[R z] 


The new result overwrites the old. 
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It IS often good practice to process the data and form [R z] before 
including the effects of a priori. When this is done one can analyze 
the effect of different a priori, [R^ z^] without reprocessing the data. 

If a priori is given in the form of an information matrix. A, 

(as for example would he the case if the problem is being initialized 

A 

with data processed using normal equation data accumulation ) then one 
can obtain R^ from A using INF2R; 


A 


INF2R 


If there were a normal equation estimate z = A b, then z — R z. 

^ o o 

II . 3 Batch Sequential Data Processing 

Prime reasons for batch sequential data processing are that many 
problems are too large to fit in core, are too expensive in terms of core 
cost, and for certain problems it is desirable to be able to incorporate 
new data as it becomes available. Subroutines TTH and UDMES are specially 
designed for this kind of problem. Both of these subroutines overwrite 
the a priori with the result which then acts as a priori for the next 
batch of data. If the data is stored on a file or tape as z^, A^, 
then the sequential process can be represented as follows: 

SRIF Processing" * 

a) Initialize [R z] with a-priori values or zero 

b) Read the next [A z] from the file 


~ — — T ^ T T 

”i.e., solving Ax = b-v with normal equations, A Ax^ -Ab; A-AA 

is the information matrix. 

A A 

The acronym SRIF represents Square Root Information Filter. The SRIF is 
discussed at length in reference [3]. 


7 



77-26 


c) 


[R z]l 
[A z]) 


THH 


[R z] 


d) If there is more data go back to b) 

e) Compute estimates and/or covariances using UTINV and RI2COV 
(as in example II. 1) 


a") 

bO 


H-D** Processing 

Initialize [U-D x] with a priori U-D information and estimate 
Read the next [A z] scalar measurement from the file 


[U-D x]| 
[A z3 ) 


[U-D x] 


d"*) If there is more data go back to b"') 

e"') Compute standard deviations or covariances using UD2SIG or 


UD2C0V. 


Note that subroutine THH is best (most efficiently) used with 
data batches of substantial size (say 5 or more) and that UDMES processes 
measurement vectors one component at a time. If the dimension of the 
state is small the cost of using either method is generally negligible. 
The UDMES subroutine is best used in problems where estimates are 
wanted with great frequency or where one wishes to monitor the effects 
of each update. In a given application one might choose to process 
data in batches for awhile and during critical periods it may be 


The new result overwrites the old. 

A* 

U-D processing is a mmerically stable algorithmic formulation of the Kalman 
filter measurement update algorithm, cf reference [3], The estimate error 
covariance is used in its UDU^ factored form, where U is unit upper triangular 
and D is diagonal. 


8 



77-26 


desirable to monitor the updating process on a point by point basis. 

In cases such as this, one may use RUDR to convert a SRIF array to TJ-D 
form or vice-versa. 

Remarks : Another case where an R to U-D conversion can be useful occurs 

in large order problems (with say 100 or more parameters) where after 
data has been SRIF processed one wants to examine estimate and/or 
covariance sensitivity to the a priori variances of only a few of the 
variables. Here it may be more convenient to update using the UDMES 
subroutine. 

11.4 Reduced State Estimates and/or Covariances From a SRIF Array 

Suppose, for example, that data has been processed and that we have a 
A A 

triangular SRIF array [R z] corresponding to the 14 parameter names, a^, a^, 

a , X, V, z, V , V , V , GM, CU41, L041, CU43, L043 (constant spacecraft 
y’ » J ’ X y z 

accelerations, position and velocity, target body gravitational constant, 
and spin axis and longitude station location errors). 

Let us ask first what would the computed error covariance be of 
a model containing only the first 10 variables, i.e., by ignoring the 
effect of the station location errors. One would apply UTINV and RI2COV 
just as in example II. 1, except here we would use N (the dimension of 
the filter ) = 10, instead of N=14. 

Next, suppose that we want the solution and associated covariance 
of the model without the 3 acceleration errors. One ESP solution is to 
use 
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r; ''t R2A rAl 
[R z] ^[A] 


NAME ORDER OF A 


X, y, z, v^, Vy, v^, 

GM, CU41-, L041, CU43, L043, 


RHS , a , a , a 
r X y 


Remark : One could also have used subroutine COMBO, with the desired 

namelist as simply a , a , a . This would achieve the same A matrix 

r’ X y 


form. 


• [A]-^SS^[R] 


Remark : R here can replace the original R and z. 

® [R]-^^[R-"x^^^]^i^[COVx^^^] 

Remarks : Here, use only N=ll, i.e., 11 variables and the RHS. 

the 11 state estimate based on a model that does not contain acceleration 

errors a , a , or a . 

r X y 

Note how triangular! zing the rearranged R matrix produces the 
desired lower dimensional SRIF array; and this is the same result one 
would obtain if the original data had been fxt using the 11 state model. 

As the last subcase of this example suppose that one is only 
interested in the SRIF array corresponding to the position and velocity 
variables. The difference between this example and the one above is 
that here we want to include the effects due to the other variables. 


z is often gxven the label RHS (right hand side) 
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One might want this sub-array to combine with a position-velocity SRIF 


array obtained from, say, optical data. One method to use would be, 

R2RA 


• [R z] 
INPUT NAMES: 




OUTPUT NAMES: 


a , a , a , X, y, z, V , V , V , GM x, y, z, v , v , v , GM 

r x y xyz xyz 

CU41, L041, CU43, L043, RHS GU41, L041, CU43, L043, RHS 

Remark ; The lower triangle starting with x is copied into R^. 

T?? A 

• [R^ Z^] ^[A, Z^l (Reordering) 

NAMES: GM, CU41, L041, CU43, L043, 

X, y, z, V , V , V , RHS 
xyz 

THH ^ ^ 

• [A, z^] ^[R^ z^] (Triangularizing) 

^ ■R9‘RA 

® [R. z ] ►[R z ] (Shifting array) 

A. A X 3IC 

NAMES: x, y, z, v , v , v , RHS 

X y z 

Remark: The lower right triangle starting with x is copied into R . 

We note that one could have elected to use COMBO in place of the first 
R2RA usage and R2A; this would have involved slightly more storage, but 
a lesser number of inputs. The sequence of operations is in this case. 


COMBO r. n 
[R z] -[A z] 


ORIGINAL NAMES 


DESIRED NAMES: x, y, z, v , v , v , RHS 

xyz 


Remark ; By using COMBO the columns of IR z] are ordered corresponding to 


the names a , a , a , GM, CU41, L041, CU43, and L043, followed by the 
r X y 

desired names list. 


11 



77-26 


TTJT.J ^ ^ 

• [A z]-!^[R z] 

Remark : The [R z] array that is output from this procedure is 

equivalent but different from the [R z] array that we began with. 


[R z]- 


-[R z ] 

K X 


Remark : As before, the lower right triangle starting with x is copied 

into R . 

X 


To delete the last k parameters from a SRIF array, it is not 
necessary to use subroutines R2A and THH. The first N - k = N columns 
of the array already correspond to a square root information matrix of 
the reduced system. If estimates are involved one can simply move the 
z column left using: 

R (N*(N + l)/2 + i) = R(N*(N +1) /2 + i) , i = l,...,k. 

Remark : We mention in passing that if one is only interested in estimates 

and/or covariances corresponding to the last k parameters then one can use 

R2RA to transform the lower right triangle of the SRIF array to an upper 

left triangle after which UTINV and RI2C0V can be applied. 

II. 5 Sensitivity, Perturbation, Computed Covariance and Consider 
Covariance Matrix Computation 

Suppose that one is given a SRIF array 



(II. 5a) 
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in X'jhich the variables are to be considered. (One can, of course, using 
subroutines R2A and THE reorder and retriangularize an arbitrarily arranged 
SRIR array so that a given set of variables fall at the end.) For various 
reasons one -may choose to ignore the y variables in the equation 


Rx+R y=z -u , V eN(0,l) 
X xy X X X 


(II. 5b) 


and take as the estimate x = R ^ z . It then follows that 

c XX 


-1 -1 

x-x=-R R y-R V , 
c X xy XX 


(II. 5c) 


and from this one obtains 


Sen = 


3 (x-x ) 



-R~^ R 
X xy 


(II. 5d) 


(sensitivity of the estimate error to the unmodeled y parameters) 

Pert = Sen Diag (o^ (1) , . . . (N^) ) (II. 5e) 

where (1) , . . . ,a^(N^) are a priori y parameter uncertainties. 

(The perturbations are a measure of how much the estimate error could be 
expected to change due to the unmodeled y parameters.) 


P =R^R^+SenP Sen^ 
con XX y 

= P^ + (Pert) (Pert) ^ if P^ is diagonal* 


(II. 5f) 


where P is the estimate error covariance of the reduced model. 


c 

An easy way to compute P^, Pert 
R2RA to place the y variable a priori 


and P 

con 



is as follows: Use subroutine 

Yo]** into the lower right 


* k 

Pert = Sen P 

y 

The a priori estimate y^ of consider parameters is generally zero. 
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corner of (II. 5a), replacing R and z , i.e.. 



Now apply subroutine UTIROW to this system (with a -1 set in the lower right 
comer* ) 



Note that the lower portion of the matrix is left unaltered, i.e., the purpose 

of UTIROW is to invert a triangular matrix, given that the lower rows have 

already been inverted. From this array one can, using subroutine RI2C0V, 

get both P and P 

c con 

c 

j-j^-lj _R12C0V^ j-p^] computed covariance 



Pert] 


RI2C0V 


[P ] 

con-^ 


consider covariance 


Suppose now that one is dealing with a U-D factored Kalman filter for- 
mulation. In this case estimate error sensitivities can be sequentially 

— . 

k 

To have estimates from the triangular inversion routines one sets a -1 in the 
last column (below the right hand side) . 

Strictly speaking this is not what we call the perturbation unless Ry(0) is 
diagonal . 
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calculated as each scalar measurement 


T T 

(2 = a X + a y + v) 
X y 


IS processed. 


Sen = Sen, , - K. (a^ Sen + a'^) 

J J-1 J X 3-1 y" 


where Sen^ ^ is the sensitivity prior to processing this (j-th) measurement, 

and K. is the Kalman gain vector. In this formulation one computes P 

j con 

in a manner analogous to that described in section II. 7; 

Let U. = U, , D- = D. (filter U-D factors) 

1 J 1' 1 

[s , ] “ Sen. (estimate error sensitivities) 
then compute 




AGTRN 


U, , _ -D, - 
k+1 k+1 


k - 1, . . , , n 


For the final U-L we have 


= U , = D 

1+1 n +1 1+1 n +1 

y y 


If P (0) = U D U j instead of P (0) = Diag (cr. , . . . , a ), then in the 

y y y y y ^ i’ ’ 

U-D recursion one should replace the Sen. columns by those of Sen.U_, and 

J J J 

0^ should be replaced by the corresponding diagonal elements of D • 

3 ^ 

II . 6 Combining Various Data Sets 

In this example we collect several related problems involving data sets 
with different parameter lists. 

Suppose that the parameter namelist of the current data does not 
correspond to that of the a priori SEIF array. If the new data involves 
a permutation or a subset of the SEIF namelist then an application of 
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subroutxne PEBMUT wxll create the desired data rearrangement. If the data 
involves parameters not present in the SEII’ namelist then one could use 
sub:^outxne R2A to modify the SRIF array to include the nexj names and then 
if necessary use PEEMUT on the data, to rearrange it compatibly. 

Suppose now that two data sets are to be combined and that each 
contains parameters peculiar to it (and of course there are common para- 
meters) , For example let data set 1 contain names ABC and data set 2 
contain names DEB, One could handle such a problem by noting that the list 
ABODE contains both name lists. Thus one could use subroutine PERMUT 
on each data set comparing it to the master list, ABODE, and then the 
results could be combined using subroutine TIIH An alternative automated 
method for handling this problem is to use sub routine COMBO with data 
set 1 (assuming it is in triangular form) and lamelist 2 . The result 
would be data set 1 in double subscripted form and arranged to the name- 
list AGDEB (names A and C are peculiar to dats set 1 and are put first). 
Having determined the namelist one could applj subroutine PERMUT to data 
set 2 and give it a compatible namelist ordering. 

The process of increasing the namelist stze to accommodate new 
variables can lead to problems with excessively long namelists, i.e., 
with high dimension. If it is known that a certain set of variables 
will not occur in future data sets then these variables can be eliminated 
and the problem dimension reduced. To eliminate a vector y from a SRIF 
array, first use subroutine E2A to put the y names fi-Tst in the namelisc; 
then use subroutine THH to retriangularize a-id finally use subroutine EZRA 
to put the y independent subarray in position for further use; viz. 


16 



77-26 




R R z 

y yx y 


R2RA 


[R z ] 

XX 


0 R 


X 


The rows [R^ R^^ z^] can be used to recover a y estimate (and its covariance) 
when an estimate for x (and xts covariance) are determined. (See example 
II. 4) . 


Still another application related to the combxning of data sets involves 
the combining of SRIF triangular data arrays. One might encounter such prob- 
lems when combining data from different space mxssions (that involve common 
parameters) or one might choose to process data of each type* or tracking 
station separately and then combxne the resultxng SRIF arrays. Triangular 
arrays can be combxne d using subroutine TTHH, assuming that subroutines 
R2A, THH and R2RA have been used previously to formulate a common parameter 
set for each of the sub problems. 


II. 7 Batch Sequential Whxte Noise 

It is not uncommon to have a problem where each data set contains a 
set of parameters that ^pply only to that set and not to any other, viz. 
the data is of the form 


A,x+B,y. - z . -V. 
3 3 3 3 3 


j = 


,N 


where there is generally a priori Information on the vector y^ variables. 
Rather than form a concatenated state vector composed of x, y^,...,y^ 
which might create a problem involving exhorbitant amounts of storage and 
computatxon we solve the problem as follows.' Apply subroutine THH to 

Aj^ wxth the corresponding R initialized with the y^ a priori. The 

res-ulting SRIF array is of the form 


viz. range, doppler, optical, etc. 
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N 


R 


R 


^2^ 


X, 


X, 


Copy the top N rows if one will later want an estimate or covariance of 

the y parameters. Apply subroutine TZERO to zero the top N rows and 

* 

using subroutine R2RA set in the a priori . This SRIF array is now 

ready to be combined with the second set of data [B^ and the procedure 

repeated. 

A somewhat analogous situation is represented by the class of problems 
that involve noisy model variations, i.e., the state at step j+l satisfies 


X 


i+1 


X. + G. w. 
1 J 2 


where matrix G. is defined so that w is Independent of x and w.eH(0,Q.). 

2 2 J J J 

Models of this type are used to reflect that the problem at hand is not 
truly one of parameter estimation, and that sone (or all) of the components 
vary in >a random (or at least unknown) manner that is statistically 
bounded. To solve this problem in a SRIF foritulation suppose that a priori 
for x^ and w^ are written in data equation foi'm (cf ref. [3]), 


R. X = z - V 
3 3 3 3 


v.eN(0,I ) 
1 n 

w 



v^^"^£N(0,l) 


where Qf is a Cholesky factor 


of Q. that is obtainable from C0V2RI. 


Combining 


these two equations with the one for x^_^^ 


glares 


\ 

In this example it is assumed that all of the y. variables have the same 
dimension. This assumption, though not essentill, simplifies our description 
of the procedure. 


18 



77-26 


I 

n 

w 

0 


’a 

w . 
J 


0 

-R G.Q^ 
111 

t— 


""j+l 


2 . 
J 




THH, i.e.. 


This is the equation to be trxangularized with subroutine 


Dim w Dim x 1 


Dim w { 

I 

n 

w 

0 

0 

THH ^ 

rR<"> 

1 

1 

W 

z 

1 

Dim X { 

% 

-R G Q r 
_ 111 

R 

1 

z . 


0 

^1+1 

^1+1_ 


I£ the problem is arranged so that is diagonal one can reduce storage and 
computation. The form of this algorithm is designed to allow the use of 


singular matrices. 

When the a priori for x and Q, are given in U— D factored form, 
one can obtain the U-D factors for follows: 

Let Q = (use C0V2TJD if necessary) 

j 

Set G = G = [g.,..., = Diag(d^ d^ ) 

^ w 

Apply subroutine AGTRN n^ times, with , D^ = D^ 


w 


(M) 


k+1 






k l,...,n 


w 


Then U. = U , D. = D 
j+1 n^ 1+1 
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Certain filtering problems involve dynamic models of the form 


’^j+l 


$. X. + G. w 
J 3 3 3 


Given an estimate for x., x., the predicted estimate for x. denoted 

3 3 3+1 

x^^^ is simply 



A 

X. 

3 


The U-D factors of the estimate error corresponding to the estimate 
can be obtained using the weighted Gram-Schmidt triangularization 
subroutine 

g] , Dlag - D.^P 


II. 8 Miscellaneous Uses of the Various ESP Subroutines 

In certain parameter analyses we may want to reprocess a set of data 

suppressing different subsets of variables. In this case the oiiginal 

/' 

data should be left unaltered and subroutxne 7 2A1 used to copy A into A^, 
which then can be modified as dictated by rue analysis. 

Covariance analyses sometimes are initxa'-ized using a covariance 
matrix from a different pioblem for a differed phase of the same problem) 
In such cases it may be necessary to permute, delete or insert rows and 
columns into the covariance matrix; and that can be achieved using sub- 
routine C2C. 

If a priori for the problem at hand is given as a covariance matrix 
then one can compute the corresponding SRIF or U-D initialization using 


In statistical notation that is commonly used, one writes 

xfl+ljl) = x(j jj) 
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subroutines C0V2E.I or C0V2UD. Of course, if the covariance is diagonal 
the appropriate R and U-D factors can be obtained more simply. To 
convert a priori given in the form of an information matrix to a corres- 
ponding SRIF matrix one applies subroutine INF2R. To display covariance 
results corresponding to the SRIF or U-D filter one can use subroutines 
UTINV, R12C0V and UD2C0V. The vector stored covariance results are 
displayed in a triangular format using subroutine TRIMAT. 

As ide : After careful consideration it was decided that subroutines to 

multiply matrices would not be included in our ESP. Our reasons are 
that parameter estimation does not, in the main, involve matrix 
multiplication; and when such products occur they generally involve 
matrices with special structures (viz. rectangle x triangle, triangle x 
rectangle, diagonal x triangle, etc). To see that these computations 
are not lengthy or complicated we illustrate how to compute z = Ex 
where R is a triangular vector stored matrix and x is an N vector, 

11=0 

DO 2 1=1, W 

SUM=0. 

II=II+I @11= (I, I) 

IK=II 

DO 1 K=I,N 

SUM=SUM+R(IK)*x(K) @IK=(I,K) 

1 IK=IK+K 

2 z(I)=SUM @z can overwrite x if desired. 
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Note that the II and IK incremental recursions are used to circumvent 
the W(N+l)/2 calculations of IK=K(K-1) /2+I. 

A later more encyclopedxc subroutine directory may include the 
various matrix products that occur in linear algebra applications. 

End of Aside 
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III. SUBROUTINE DIRECTORY SUMMARY 

1. AGTRN - (Agee— Turner) 

Computes updated U-D factors corresponding to a rank 1 matrix 

modificatxon; i.e., given U-D, a scalar c, and vector v, U and D are 

T T T 

computed so that UDU =UDU +cvv. Bothc and v are destroyed during 

the computation, and the resultant (vector stored) U-D array replaces 

the original one. Uses for this routine include (a) adding process 

noise effects to a U-D factored Kalman filter; (b) computing consider 

covariances (cf Section II. 5); (c) computing "actual" covariance 

factors resulting from the use of suboptimal Kalman filter gains; and 

(d) adding measurements to a U-D factored information matrix. 

2. A2A1 - (A to Al) 

Reorders the columns of a rectangular matrix A, storing the 
result in matrix Al, Columns can be deleted and new columns added. 

Zero columns are inserted which correspond to new column name entries. 
Matrices A and Al cannot share common storage. 

Example III.l 


a 

B 

C 

B 

F 

G 

C 

H 

1 

5 

9 


5 

0 

0 

9 

0 

2 

6 

10 

A2A1 

6 

0 

0 

10 

0 

3 

7 

11 


7 

0 

0 

11 

0 

-A 

8 

12 


8 

0 

0 

12 

0 


A Al 

The new namelist (BFGCH) contains F, G and H as new columns and deletes 
the column corresponding to name a. 
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Example III. 2 

Suppose one is given an observation data file" with regression 

coefficients corresponding to a state vector with components say, 

X, y, z, V , V , V and station location errors. Suppose further, 
X y z ’ 

* A A 


that the vector being estimated has components a^, a^ , 

X, y, z, V , V , V , GM and station location errors. A2A1 can be used 
X y z 

to reorder the matrix of regression coefficients to correspond to the 
state being estimated. Zero coefficients are set in-place for the 
accelerations and GM which are not present in the original file. 


3. COMBO - (combine R and A namelists) 

The upper triangular vector stored matrix R has its columns 
permuted and is copied into matrix A. The names associated with R 
are to be combined with a second namelist. 

The namelist for A is arranged so that R names not contained in 
the second list appear first (left most) . These are then followed by 
the second list. Names in the second list that do not appear in the 
R namelist have columns of zeros associated with them. 

Example HI. 3 

NAM2 list 


a 

B 

C 

D 

C 

B 

E 

a 

F 

D 

“ 1 

2 

4 

7“ 


"4 

2 

0 

1 

0 

7 " 

0 

3 

5 

8 



3 

0 

0 

0 

8 

0 

0 

6 

9 



0 

0 

0 

0 

9 

0 

0 

0 

10 



0 

0 

0 

0 

10 


R— Vector stored A-Double subscripted 


A 

in track and cross track accelerations 
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A principal application of this subroutine is to the problem of 
combining equation sets containing different variables, and automating 
the process of combxning name lists. 

4. C0V2R1 - (Covariance to E. inverse) 

An input positive semi-definite vector stored matrix P is replaced by 

T 

xts upper triangular vector stored Cholesky factor U, P = UU . The name RI 

-1 

is used because when the input covariance is positive definite, U = R 

5. C0V2UD - (Covariance to U-D factors) 

An input positive semi-definite vector stored matrix P is re- 

T 

placed by its upper triangular vector stored TJ-D factors. P = UDU . 

6. C^ - (C to C) 

Reorders the rows and columns of a square (double subscripted) 
matrix C and stores the result back in C, Rows and columns of zeros 
are added when new column entries are added. 

Example III. 4 



A 

B 

r 



r 

p 

B 

Q 

A 

1 

4 

7 


r 

9 

0 

6 

0 

B 

2 

5 

8 

C2C 

p 

0 

0 

0 

0 

r 

3 

6 

9 


B 

8 

0 

5 

0 






Q 

0 

0 

0 

0 


Names P and Q have been added and name A deleted. An Important appli 

cation of this subroutine is to the rearranging of covariance matrices. 

7. INF2R - (Information matrix to R) 

Replaces a vector stored positive semi-definite information matrix 

T T 

A by its lower triangular Cholesky factor R ; A = R R. The upper tri- 
angular matrix R is in the form utilized by the SRIP algorithms. The 
algorithm is designed to handle singular matrices because it is a 
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common practice to omit a priori information on parameters that are 
either poorly known or which will be well determined by the data. 

8. PERMUT 

Reorders the columns of matrix A, storing the result back in A. 
This routine differs from A2A1 principally in that here the result 
overwrites A. PERMUT is especially useful in applications where 
storage is at a premium or where the problem is of a recursive nature. 

9. RINCON - (R inverse with condition number bound, CNB) 

Computes the inverse of an upper triangular vector stored matrix R 
using subroutine UTINV. A Frobenius bound (CNB) for the condition 
number of R is computed too. This bound acts as both an upper and a 
lower bound, because CNB/N < condition number < CNB. When this bound is 
within several orders of magnitude of the machine accuracy the computed 
inverse is not to be trusted, (viz if CNB ^ 10 ^ on an 18 decimal digit 
machine R is ill-conditioned) . 

RI2C0V - (RI to covariance) 

This subroutine computes sigmas (standard deviations) and/or the 
covariance of a vector stored upper triangular square root covariance 
matrix, RINV (SRIF inverse). The result, stored in COVOUT (covariance 
output) is also vector stored, COVOUT can overwrite RINV. 

11. R^ - (R to A) 

The columns of a vector stored upper triangular matrix R are per- 
muted and variables are added and/or deleted. The result is stored in 
the double subscripted matrix A. In other respects the subroutine is 


like A2A1. 
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Example III. 5 


a 

B 

C 

D 

E 


E 

F 

c 

B 

" 2 

4 

8 

14 

22 " 

* 

22 

0 

8 

4 “ 

0 

6 

10 

16 

24 


24 

0 

10 

6 

0 

0 

12 

18 

26 

R2A 

26 

0 

12 

0 

0 

0 

0 

20 

28 


28 

0 

0 

0 

0 

0 

0 

0 

30 


30 

0 

0 

0 


R A 

R IS vector stored as R = (2,4,6,8,10,12,14,16,18,20,22,24,26,28,30) 
with namelist (a,B,C,D,E) associated with it. Names a and D are 
not included in matrix A, and a column of zeros corresponding to name 
F is added. 

One trivial, but perhaps useful, application is to convert a 

* 

vector stored matrix to a double subscripted form. R2A is used most 

often when one wants to rearrange the columns of a SRIF array so that 

» 

reduced order estimates, sensitivites, etc. can be obtained; or so that 
data sets containing different parameters can be combined. 

12. R2RA - (Triangular block of R to triangular block of RA) 

A triangular portion of the vector stored upper triangular matrix R 
IS put into a triangular portion of the vector stored matrix RA. The 
names corresponding to the relocated block are also moved. R can 
coincide with RA. 

A 

see also the aside in the introduction 
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Examples III. 6 


Q Z Q Z 



Note that an upper left triangular submatrix can slide to any lower 
position along the diagonal, but that a submatrix moving up must go 
to the upper leftmost corner. Upper shifting is used when one' is 
interested in that subsystem; and the lower shifting is used, for 

example, when inserting a priori information for consider analyses. 

I 

13. RUDR - (SRIF R converted to U-D form or vice versa) 

A vector stored SRIP array is replaced by a vector stored U— D 
form or conversely. A point to be noted is that when data is involved 
the right side of the SRIF data equation transforms to the estimate 
in the U-D array. 
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14. THH - (Triangular Householder data packing) 

An upper triangular vector stored matrix R is combined with a rec- 
tangular doubly subscripted matrix A by means of Householder orthogonal 
transformations. The result overwrites R, and A is destroyed in the process. 



15. TRIMAT — (Triangular matrix print) 

Prints a vector stored upper triangular matrix, using a matrix 
format. 

Example III. 7 

R(10) = (2,4,6,8,10,12,14,16,18,20) with assocxated namelist 
(A,B,C,D) is printed as 

A B C D 

A 2 4 8 14 

B 6 10 16 

C 12 18 

D 20 

(The numbers are printed to 8 significant floating point 
digits) . 

To appreciate the importance of this subroutine compare the vector 
R(10) with the double subscript representation. 

16. TTHH - (Two triangular arrays are combined using Householder 

orthogonal transformations) 

This subroutine combines two single subscripted upper triangular 
SRIE arrays, R and RA using Householder orthogonal transformations. 

The result overwrites R. 

* 

The elements are not explicitly set to zero. 
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17. TZERO - (Zero a horizontal segment of a vector stored upper 
triangular matrix) 

Upper triangular vector stored matrix R has its rows between ISTART 
and IFINAL set to zero. 

Example III. 8 

To zero row 2 and 3 of R(15), in the example of subroutine 11- 
R(15) = (2,4,6,8,10,12,14,16,18,20,22,24,26,28,30) 

R(15) = (2,4,0,8,0,0,14,0,0,20,22,0,0,28,30) 


2 

4 

8 

14 

22 


2 

4 

8 ' 

14 

22 

0 

6 

10 

16 

24 


0 

0 

0 

0 

0 

0 

0 

12 

18 

26 

TZERO^ 

0 

0 

0 

0 

0 

0 

0 

0 

20 

28 


0 

0 

0 

20 

28 

0 

0 

0 

0 

30 


0 

0 

0 

0 

30 


R-vector stored R-vector stored 


The elements are not explicitly set to zero. 
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18. UDMES - (TJ-D measurement update) 

Given the U-D factors of the a priori estimate error covariance 
and the measurement, z = Ax + v this routine computes the updated 
estimate and U-D covariance factors, the predicted residual, the 
predicted residual variance, and the normalized Kalman gain. This 
is Bierman’s U-D measurement update algorithm. 

19. UD2C0V - (U-D factors to covariance) 

The input vector stored U-D matrix (diagonal D elements are 

stored as the diagonal entries of U) is replaced by the covariance P, 

T 

also vector stored. P = UDU . P can overwrite U to economize on storage. 

20. UD2SIG - CU“D factors to sigmas) 

Standard deviations corresponding to the diagonal elements of the 

covariance are computed from the U-D factors. This subroutine, a restricted 
version of UD2C0V can print out the resulting sigmas and a title. The 
input U-D matrix is unaltered. 

21. UTINV - (Upper triangular matrix Inversion) 

An upper triangular vector stored matrix RIN(R in) is inverted and 
the result, vector stored, is put in R0UT(R out). ROUT can overwrite 
RIN to economize on storage. If a right hand side is included and the 
bottommost tip of RIN has a -1 set in then ROUT will have the solution 
in the place of the right hand side. 
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22. UTIROW - (Upper triangular Inversion, inverting only the upper rows) 


INPUT OUTPUT 



R 

X 

1 

UTIROW 

R ^ 

X 

-r“^r R ^ 

X xy y 

n 

0 

R“1 

0 

R-^ 

y 


y 



y 


An input vector stored R matrix with its lower left triangle assumed to 
have been already inverted is used to construct the upper rows of the 
matrix inverse of the result. The result, vector stored, can overwrite 
the input to economize on storage. 

If the columns comprising R^ represent consider terms then taking 
R_^ as the identity gives the sensitivity on the upper right portion of 

y ^ 

-1 

the result. If R = Diag(a , . . . ,a ) then the upper right portion of 
y y 

the result represents the perturbation. Note that if z (the right hand 
side of the data equation) is included in R^ then taking the corres- 
ponding R^^ diagonal as -1 results in the filter estimate appearing 
as the corresponding column of the output array, l^hen n^ is zero this 
subroutine is equivalent to UTINV, 

23. WGS - CWeighted Gram Schmidt matrix triangularization) 

An input rectangular (possibly square) matrix W and a diagonal 
weight matrix, D^, are transformed to (U-D) form; i.e,, 

T T 

S D W = UDU 
w 

where U is unit upper triangular and D is diagonal. The weigh,ts are 
assumed nonnegative, and this characteristic is inherited by the 
resulting D, 
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IV. SUBROUTINE DIRECTORY USER DESCRIPTION 

1. AGTRN (Agee-Tumer U-D rank one modification) 

Purpose 

T T 

Xo compute the (updated) U-D factors of UDU + CW . 
CALL AGTBN (UIN,UOUT,N,C,V) 


Argument Definitions 

UIN(N*(N+1) /2) Input vector stored positive semi- 

definite U-D array (with the D entries 
stored on the diagonal of U) 

U0DT(N*(N+l)/2) Output vector stored result 

U0UT=UIN is allowed 


N 

C 

V(N) 


Matrix dimension 

Input scalar, destroyed by the algorithm 
Input vector, destroyed by the algorithm 


Remarks and Restrictions 

If C negative is used the algorithm is numerically unstable, 
and the result may be numerically unreliable. Singular U matrices 
are allowed, and these can result in singular output U matrices. 
Functional Description 

This rank one modification is based on a result published by 
Agee and Turner (1972), PJhite Sands Missile Range Tech. Report 
No. 38. See also Ref. [3] where the algorithm is derived using 
geometric arguments. 
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2. A2A1 (A to Al) 

Purpose 

To rearrange the columns of a namelist indexed matrix to 
conform to a desired namelist. 


CALL A2A1(A,IA,IR,LA,MAMA,A1,IA1,LA1,NAMA1) 


Argument Definitions 
A(IR,LA) 
lA 
IR 

LA 


Input rectangular matrix 

Row dimension of A, IA..GE.IR 

Number of rows of A that are to be 
arranged 

Number of columns in A; this also 
represents the number of parameter 
names associated with A 


NAMA(LA) • 

A1(IR,LA1) 

lAl 


Parameter names associated with A 

Output rectangular matrix 

Row dimension of Al, lAl.GE.lR 


LAI Number of columns in Al ; this also 

represents the number of parameter 
names associated with Al 


NAMAl(LAl) 


Input list of parameter names to be 
associated with the output matrix Al 


Remarks and Restrictions 


Al cannot overwrite A. This subroutine can be used to add 

on columns corresponding to new names and/or to delete variables 

from an array . 

Functional Description 

The columns of A are copied into Al in an order corresponding 
to the NAMAl parameter namelist. Columns of zeros are inserted 
in those Al columns which do not correspond to names in the input 
parameter namelist NAMA. 
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3. COMBO (Combine parameter namelists) 
Purpose 


To rearrange a vector stored triangular matrix and store 
the result in matrix A. The difference between this subroutine 


and R2A is that there the namelist for A is input 5 here it is 
determined by combining the list for R with a list of desired names. 

CALL COMBO (R,L1,NAM1,L2,NAM2,A,IA,LA,NAMA) 

Argument Definitions 


R(Ll*(Ll+l)/2) 

LI 

NAMl(Ll) 

L2 

NAM2 (L2) 

A (LI, LA) 


Input vector stored upper triangular matrix 

No. of parameters in R (and in NAMl) 

Names associated with R 

No . of parameters in NAM2 

Parameter names that are to be combined 
with R (NAMl list); these names may or 
may not be in NAMl 

Output array containing the rearranged 
R matrix Ll.LE.IA * 


Ik 


Row dimension of A 


LA 


NAMA(LA) 


No. of parameter names in NAMA, and the 
column dimension of A. LA==Ll + L2 - 
No. names common to NAMl and NAM2; LA 
is computed and output 

Parameter names associated with the out- 
put A matrix ; consists of names in NAMl 
not in NAM2 followed by NAM2 


Remarks and Restrictions 


The column dimension of A is a result of this subroutine. 

To avoid having A overwrite neighboring arrays one can bound the 
column dimension of A by L1 + L2. 
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Functional Description 

First the NAMl and NAM2 lists are compared and the names 
appearing in NAMl only have their corresponding R column entries 
stored in A (e.g. if NAM1(2) and NAM1(6) are the only names not 
appearing in the NAM2 list then columns 2 and 6 of R are copied 
into columns 1 and 2 of A) . The remaining columns of A are 
labeled with NAM2. The A namelist is recorded in NAMA, The 
NAMl list is compared with NAM2 and matching names have their R 
column entries copied into the appropriate columns of A. NAM2 
entries not appearing in NAMl have columns of zero placed in A. 
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4. C0V2RI (Covariance to Cholesky Square Root, RI) 

Purpose 

To construct the upper triangular Cholesky factors of a 
positive semi-definite matrix. Both the input covariance and 
the output Cholesky factor (square root) are vector stored. 

The output overwrites the input. Covariance (input) = U*U**T 
(output U = Einverse) . 

CALL C0V2RI(U,N) 

Argument Definitions 

U(N*(N+l)/2) Contains the input vector stored 

covariance matrix (assumed positive 
definite) and on output it contains 
the upper triangular square root factor 

N Dimension of the matrices involved 

Remarks and Restrictions 

No check is made that the input matrix is positive semi- 
definite. Singular factors (with zero columns) are obtained If 
the input is (a) in fact singular, (b) ill-conditioned, or (c) in 
fact indefinite; and the latter two situations are cause for alarm. 
Case (c) and possibly (b) can be identified by using RI2C0V to 
reconstruct the input matrix. 

Fmctional Description 

An upper triangular Cholesky reduction of the input matrix 

is implemented using a geometric algorithm described in Ref. [3]. 

T 

U( input) = U( output) * U( output) 

At each step of the reduction diagonal testing is used and 
negative terms are set to zero. 
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5. C0V2UD (Covariance to UD factors) 

Purpose 

To obtain the U-D factors of a positive semi— definite matrix. 

The input vector stored matrix is overwritten by the output U-D 
factors which are also vector stored. 

~CALL C0V2UD(U~ N) 

Argument Definitions 

U(N*(N+l)/2) Contains the input vector stored covari- 

ance matrix; on output it contains the 
vector stored U-D covariance factors. 

N Matrix dimension 

Remarks and Restrictions 

No checks are made in this routine to test that the input U matrix 
is positive semi-definite. Singular results (with zero columns) are 
obtained if the input is (a) in fact singular, (b) ill-conditioned, 
or (c) in fact indefinite; and the latter two situations are cause for 
alarm. Case (c) and possibly case (b) can be identified by using UD2- 
COV to reconstruct the input matrix. Note that although indefinite 
matrices have D-D factorizations, the algorithm here applies only to 
matrices with non-negative eigenvalues . 

Functional Description 

An upper triangular U-D Cholesky factorization of the input matrix 

is implemented using a geometric algorithm described in Ref. [3]. 

T 

U(input) = U*D*U , U-D stored in U on output 

at each step of the reduction diagonal testing is used to zero negative 
terms . 
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6. C2C (C to C) 

Purpose 

To rearrange the rows and columns of C, from NAMl order to NAM2 
order. Zero rows and columns are associated with output defined names 
that are not contained in NAMl. 

CALL C2C(C,IC,L1,NAM1,L2,NAM2) 

Argument Definitions 

C(L1,L1) 

IC 

LI 

NAMl(L) 

L2 

NAM2 (L2) 

Remarks and Restrictions 

The NAM2 list need not contain all the original NAMl names and 
LI can be .GE. or .LE. L2. The NAMl list is used for scratch and 
appears permuted on output. If L2.GT.Ll the user must be sure that 
NAMl has L2 entries available for scratch purposes. 

Functional Description 

The rows and columns of C and NAMl are permuted pairwise to get 
the names common to NAMl and NAM2 to coalesce. Then the remaining rows 
and columns of C(L2,L2) are set to zero. 


Input matrix 

Row dimension of C 
IC.GE.L = MAX(L1,L2) 

No. of parameter names associated with 
the input C 

Parameter names associated with C on input . 
(Only the first LI entries apply to the 
input C) 

No. of parameter names associated with the 
output C 

Parameter names associated with the output C 
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7 . INF2R (Information matrix to R) 

Purpose 

To compute a lower triangular Cholesky factorization of the 
input positive semi-definite matrix. The result transposed, is 
vector stored; this is the form of an upper triangular SRIF matrix. 

CALL INF2R(P,N) 

Argmnent Definitions 

P(N*(N+l)/2) Input vector stored positive semi- 

definite (information) matrix; on output 
it represents the transposed lower 
triangular Cholesky factor (i.e. the SRIF 
R matrix) 

N Matrix dimension 

Remarks and Restrictions 

No checks are made on the input matrix to guard against negative 
eigenvalues of the input, or to detect ill-conditioning. Singular 
output matrices have one or more rows of zeros . 

Functional Description 

A Cholesky type lower triangular factorization of the input matrix 
is implemented using the geometric formulation described in Ref. [3]. 
Il(lnput) = [U(output) *[U(output) ] 

At each step of the factorization diagonal testing is used to zero columns 
corresponding to negative entries. The result is vector stored in the 
form of a square root information matrix as it would he used for SRIF 
analyses . 
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8. PERMUT (Permute A) 


Purpose 

To rearrange the columns of a namelist indexed matrix to conform 


to a desired namelist. The resulting matrix is to overwrite the input. 


CALL PERMUT(A,IA,IR,L1,NAM1,L2,NAM2)^ 


Argument Definitions 
A(IR,L) 
lA 
IR 

LI 

NAMl (L) 


L2 

NAM2 


Input rectangular matrix, L = max(Ll,L2) 

Row dimension of A, lA.GE.IR 

Number of rows of A that are to be 
rearranged 

Number of parameter names associated with 
the input A matrix 

Parameter names associated with A on input 
(only the first LI entries apply to the 
input A) 

Number of parameter names associated with 
the output A matrix 

Parameter names associated with the output A 


Remarks and Restrictions 


This subroutine is similar to A2A1; but because the output matrix 
in this case overwrites the input there are several differences. The 
NAMl vector is used for scratch, and on output it contains a permuta- 
tion of the input NAMl list. The user must allocate L = max(Ll,L2) 
elements of storage to NAMl. The extra entries, when L2>Ll, are 
used for scratch. 

Functional Description 

The columns of A are rearranged, a pair at a time, to match the 
NAM2 parameter namelist. The NAMl entries are permuted along with the 
columns, and this is why dim (NAMl) must be larger than LI (when L2>Ll) . 
Columns of zeroes are inserted in A which correspond to output names 
that do not appear in NAMl. 
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RINCON (R inverse with condition number bound) 

Purpose 

To compute the inverse of an upper triangular vector stored 
triangular matrix, and an estimate of its condition number. 

CALL RINCON (RIN.jN ,R0I3T., CNB) | 

Argument Definitions 
RIN(N*(N+l)/2) 

N 

R0UT(N*(N+l)/2) 

CNB 

Remarks and Restrictions 

The condition number bound, CNB serves as an estimate of the actual 
condition number. When it is large the problem is ill-conditioned. The 
matrix inversion is computed using subroutine UTINV. 

Functional Description 

The matrix inversion, a triangular back substitution, is accomplished 
via subroutine UTINV. If any diagonal element of the input R matrix is 
zero the inversion is not attempted; instead a message is printed. The 
condition number bound rs computed as follows: 

NTOT 

F.NORH R = R(J)^ 

J=1 

NTOT 

F.NORM r"^ = R-^CJ)^ 

J=1 


Input vector stored upper triangular matrix 
Matrix dimension 

Output vector stored matrix inverse 
(RIN = ROUT is permitted) 

Condition number bound. If k is the 
condition number of RIN, then 
CNB/N.LE.k.LE CNB 
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where NTOT = N*(N+-l)/2 is the number of elements in the vector stored 

r 

triangular matrix. The condition number bound, CNB, is given by 

-1 1/2 

CNB = (F.NORM R * F.NORM R ) 

F.NORM is the Frobenius norm, squared. The inequality 

CNB/N < condition number R ^ CNB 

is a simple consequence of the Frobenlus norm inequalities given in 
Lawson-Hanson "Solving Least Squares," page 234. 
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10. RI2C0V (RI Triangular to covariance) 

Purpose 

To compute the covariance matrix and/or the standard deviation of 
a vector stored upper triangular square root covariance matrix. The 
output covariance matrix, also vector stored, may overwrite the input. 
"call RI2COV(RINV,N,SIG,COVOUT,KOV) 


Input vector stored upper triangular 
covariance square root (RINV=R inverse 
is the inverse of the SRIF matrix) . 

Dimension of the RINV matrix 

Output vector of standard deviations 

Output vector stored covariance matrix 
(COVOUT = RINV is allowed) 

Compute covariance and sigmas using the 
first KOV rows of RINV 

Compute only the sigmas using the first 
KOV rows of RINV 

No covariance, but all sigmas (e.g. use 
all N rows of RINV) 


Argument Definitions 
RINV(N*(N+l)/2) 

N 

SIG(N) 

COVOUT (N*(N+1)/ 2) 

1 .GT.0 
.LT.O 
.EQ.O 


Remarks and Restrictions 

Replacing N by ]kOV| corresponds to computing the covariance of 
a lower dimensional system. 

Ftinctional Description 

COVOUT=RINV*RINV**T . 
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11. E2A (R to A) 


Purpose 

To place the upper triangular vector stored matrix R into the 
matrix A and to arrange the columns to match the desired NAMA para- 
meter list. Names in the NAMA list that do not correspond to any 
name in NAMR have zero entries in the corresponding A columns. 


CALL R2A(R,LR,NAMR,A,IA,LA,NAMA) 


Argument Definitions 
R(LR*(LR+l)/2) 
LR 

NAMR(LR) 

A(LR,LA) 

lA 

LA 

NAMA (LA) 


Input upper triangular vector stored array 

Row dimension of vector stored R 

Parameter names associated with R 

Matrix to house the rearranged R matrix 

Row dimension of A, lA.GE.LR. 

No . of parameter names associated with the 
output A matrix. 

Parameter names for the output A matrix. 


Functional Description 

The matrix A is set to zero and then the columns of R are copied 
into A. 
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12, R2EA (Permute a subportion R of a vector stored triangular matrix) 

XX 

Purpose 

To copy the upper left (lower right) portion of a vector stored 
upper triangular matrix R into the lower right (upper left) portion of 
a vector stored triangular matrix RA., 

CALL R2RA(R,NR,NAM,RA,NRA,RAMA) 


Argument Definitions 

R(NR*(NR+1) /2) Input vector stored upper triangular matrix 


NR 

NAM (NR) 

RA(NRA*(NRA+l)/2) 


Dimension of vector stored R matrix'^' 

Names associated with R. 

Output vector stored upper triangular matrix 


NRA 


NAMA(NRA) 


If NRA= 0 on input, then NAMA(l) should have 
the first name of the output namelist. In 
this case the number of names in NAMA, NRA, 
will be computed. The lower right block of 
R will be the upper left block of RA. 


If NRA = last name of the upper left block 
that is to be moved then this upper block 
is to be moved to the lower right corner 
of RA, When used in this mode NRA=NR on 
outputs 


Names associated with RA. Note that NRA 
used here denotes the output value of NRA. 


Remarks and Restrictions 


RA and NAMA can overc^rite R and NAM. The meaning of the WRA= 0 


option is clarified by the following example; 


A 

B 

C 

D 

E 

C 

D 

E 

INPUT 







- 


m 

NR = 5 

2 

4 

8 

14 

11 


12 

18 

26 

NAM = 'A’ 










NRA = 0 


6 

10 

16 

24 



20 

28 

NAMA(l) = 

R 



il2 

18 

26 




30 

OUTPUT 



1 




- 


• 

NAMA = 'C 



1 

20 

28 







1 


30_ 


RA 




R 


^see the concluding paragraph of Remarks and Restrictions 
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When NRA, = 0 and NAMA.(1) = 'C* we are asking that the lower triangular 
portion of R, beginning at the column labeled C, be moved to form the 
first (in this case 3) columns of RA.. Incidentally, RA could have 
additional columns; these columns and their names would be unaltered 
by the subroutine. 

The meaning of the other NRA option is illustrated by the following 
example; 



INPUT 
NR = 5 
NAM = 'A', 

NRA = ’C' 

R 

OUTPUT 
NRA = 5 

NAMA(3-5) = 'A','B\'C 
RA 


%en NRA = ’C’ we are asking that the upper left block of R, up to the 
column labeled C, be moved to the lower left portion of RA. and the cor- 
responding names be moved too. If RA overwrites R, as in the example, 
then the first two rows of R remain unchanged and since NAMA overwrites 
NAM, the labels of the first two columns remain unaltered. 

The remark that NRA=NR on output means, in this example, that the 
column with name C in R is moved over to column 5. If one wanted to 
slide the upper left triangle corresponding to names ABC of R to columns 
7-9 of an RA matrix (of unspecified dimension, > 9), then one should set 
NR=9 in the subroutine call. Thus NR, when used in this sliding down 
the diagonal mode, does not represent the dimension of R; but indicates 
how far the slide will be. 
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L3. RUDR (R to U-D or U-D to R) 
Pu rpose 


To transform an upper triangular vector stored SRIF array to U-D 
form or vice versa. 

CALL RUDR (RIN.N, ROUT, IS) 


Argument Definitions 

RIN (NEAR* (NBAR+1) /2) Input upper triangular vector stored SRIF 

or U-D array; NEAR = ABS(N) + 1 

R0UT(NBAR*(NBAR+1)/2) Output upper triangular vector stored 

U-D or SRIF array (RIN = ROUT is 
permitted) 


N Matrix dimension, N.GT.O represents an 

R to U-D conversion and N.LX.O represents 
a U-D to R conversion. 


IS . If IS = 0 the input array is assumed not 

to contain a right side (or an estimate) , 
and IS = 1 means an appropriate additional 
column is included. In the IS = 0 case 
the last column of RIN is ignored and 
NEAR = ABS(N) is used. 


Subroutine used: UTINV 


Functional Description 

Consider the N>0 case. RIN = R is transformed to ROUT — R inverse 

using subroutine UTINV with dimension N+IS. If IS = 1 the subroutine 

sets RIN((N+1) (N+2))/2) = -1. so that the N+lst column of ROUT will be 

-1 1/2 

the X estimate followed by —1. R = UD so that the diagonals 
are square root scaled U columns. This information is used to con- 
struct the U-D array which overwrites ROUT, 

Tf N<0 the input is assumed to be a U— arra^. This array is 
converted to ROUT=UD^ and then using UTINV, R is computed and stored in 
ROUT, If IS = 1 the U-D matrix is assumed augmented by X (estimate), 
and on output the right side term of the SRIF array is obtained. 
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14. THH (Triangular Householder Orthogonalizatlon) 


Purpose 


To compute [R z] 


such that 



T - orthogonal 


This is the key algorithm used xn the square root informatxon batch 
sequential fxlter. 


CALL THH(R,N,A,IA,M,SOS,NSTRT) 


Argument Definitions 
R(N* (N+3)/2) 


N 

A(M,N+1) 


lA 

M 

SOS 


Input upper triangular vector stored 
square root information matrix. If 
estimates are involved SOS.GE.O and R 
is augmented with the right hand side 
(stored in the last N locations of R) . 

If SOS.LT.O only the first N*(N+l)/2 
locations of R are used. The result 
of the subroutine overwrites the input R 

No. of parameters 

Input measurement matrix. The N+lst 
column is only used if SOS.GE.O, in 
which case it represents the right side 
of the equation v + AX = z. A is 
destroyed by the algorithm, but it is 
not explicitly set to zero. 

Row dimension of A 

The number of rows of A that are to be 
combined with R 

\ 

Accumulated residual sum of squares 
corresponding to the data processed 
prior to this time. On exit SOS 
represents the updated sum of squares 

of the residuals ^ I I ^ ’ 

summed over the old and new data. It 
also includes the a priori term 

I |r X “Z I t • Because SOS cannot 

be used if data, z, is not included we 
use SOS.LT.O to Indicate when data is 
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not included. 

NSTRT ' First column of the input A matrix 

that has a nonzero entry. In certain 
problems, especially those involvxng 
the inclusion of a priori statistics, 
it is known that the first NSTRT-1 
columns of A all have zero entries. 

This knowledge can be used to reduce 
computation. If nothing is known 
about A then NSTRT, LE.l gives a default 
value of 1, i.e. it is assumed that A 
may have nonzero entries in the very 
first column. 


T»er<arks and Restrictions 

It is trivial to arrange the code so that R output need not over- 
write the input R. This was not done because, in the author's opinion, 
there are too few times when one desires to have ROUT ^ RIN. 


Functional Description 

Assume for simplicity that NSTRT = 1. Then at step j, j = 1,...,N 
(or N+1 if data is present) the algorithm implicitly determines an 
elementary Householder orthogonal transformation which updates row j 
of R and all the columns of A to the right of the jth. At the 
completion of this step column j of A is in theory' zero, but it is 
not explicitly set to zero. The orthogonalization process is discussed 
at length in the books by Lawson and Hanson, [1] and Bierman [3]. 
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15 . TRIJIAT (Triangular matrix print) 

Purpose 

To display a vector stored upper triangular matrix in a two 
dimensional 8-digit triangular format. 

CALL TRIMAT (A, H, CAR, TEXT,NCHAR,NAMES) 


Argment Definitions 

A(N*N+l)/2) Vector stored upper triangular matrix 

N Dimension of A 


CAR(N) 


Parameter names (alphanumeric) associated 
with A 


TEXT(NCHAR) 

NCHAR 


NAliES 


An array of field data characters to 
be printed as a title preceding the matrix 

No. of characters (including spaces) that 
are to be printed in text( ) 

ABS (NCHAR) .LE .126 .NCHAR negative is used 
to avoid skipping to a new page to start 
printing 

A logical flag. If NAMES=.F. the CAR 
namelist is ignored and the columns 
and rows of A on output appear with 
numerical column heads 


Remarks and Restrictions 

Using NCHAR nonnegative, and starting the print at the top of a 
new page makes it easier to locate the printed result and is especially 


recommended when dealing with large dimensioned arrays. Page economy 
can, however, be achieved using the NCHAR negative option. In this case 
the print begins on the next line. 
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16 , TTHH (Two triangular matrix Householder reduction) 

Purpose 

To combine two vector stored upper t^i angular matrices , R and RA 
by applying Householder orthogonal transformations. The result over- 
writes R. 



Argument Definitions 
R(N*(N+1) /2) 


Input vector stored upper triangular 
matrix, which also houses the result 


RA(N*(N+l)/2) 


N 


Remarks and Restrictions 


Second input vector stored upper 
triangular matrix. This matrix is 
destroyed by the computation. 

Matrix dimension 

N less than zero is used to indicate 
that R and RA have right sides 
(|n|+1 columns) and have dimension 
1n|*(In|+3)/2). 


RA is theoretically zero on output, but is not set to zero. 
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17. TZERO (Triangular matrix zero) 

Purpose 

To zero out rows IS (I start) to IF(I final) of the vector stored 
upper triangular matrix R. 

CALL TZERO (R,N, IS, IF) ! 


Argument Definition 
R(N*(N+l)/2) 

N 

IS 

IF 

Functional Description 



Input vector stored upper triangular 
matrix 

Row dimension of vector stored matrix 
First row of R that is to be set to zero 
Last row of R that is to be set to zero 


IS 
IF 

R( output) 
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18. UDMES (U-D measurement update) 

Purpose 

Kalman filter measurement updating using Bierman's U=D measure- 
ment update algorithm, cf 1975 CONF. DEC. CONTROL paper. A scalar 
T 

measurement z = A x + v is processed, the covariance U-D factors and 
estimate (if included) are updated, and the Kalman gain and innovations 
variance are computed. 

CALL UDMES (U,N,R, A, G, ALPHA) 


Argument Definitions 
INPUTS 

U(N*(N+l)/2) Upper triangular vector stored xnput 

matrix. D elements are stored on the 
diagonal. The U vector corresponds to 
an a priori covariance. If state 
estimates are involved the last column 
of U contains X. In this case Dim U = 
(N+l)*(N+2)/2 and on output (U((N+1)* 
(N+2)/2) = z-A**T*X(a priori est) . 

N Dimension of the state vector 


R 

A(N) 

ALPHA 


OUTPUTS 


U 


ALPHA 


Measurement variance 

Vector of Measurement coefficients; 
if data then A(N+1) = z 

If ALPHA. LT. zero no estimates are 
computed ( and X and z need not be 
xncluded) 


Updated vector stored U-D factors. I^Then 
ALPHA (input) is nonnegative the (N+l)st 
column contains the updated estimate 
and the predicted residual. 

Innovations variance of the measurement 
residual. 


A Contains U**T*A(input) and when ALPHA 

(input) is nonnegative A(N+1) = 
z-A**T*X(a priori est) /ALPHA. 
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G(N) Vector of unweighted Kalman gains, 

K = G/ ALPHA. 

Remarks and Restrictions 

One can use this algorithm with R negative to delete a previously 
processed data point. One should, however, note that data deletion 
sometimes introduces numerical errors. 

The algorithm holds for R = 0 (a perfect measurement) but the code 
may fail (zero divides occur) if any of the ALPHA terms appearing in 
the code vanish. Changes in the code which remove the zero divide 
problems are commented in the code. 

Functional Description 

The algorithm updates the columns of U, from left to right, using 
Bierman's algorithm, cf Proc. 1975 Conf. Dec. Control, Houston, Texas, 
pp 337-346. 
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19 . UD2C0V (U-D factor to covariance) 

Purpose 

To obtain a covariance from its U-D factoriza'tion. Both matrices 

are vector stored and the output covariance can overwrite the input 

T 

U-D array. U-D and P are related via P = UDU , 

f ^ — — 

I CALL UD2C0V(UIN,N,P0UT) 


Argument Definitions 

UIN(N*(N+1) /2) Input vector stored U-D factors, with D 

entries stored on the diagonal. 

P0UT(N*(W+1) /2) Output vector stored covariance matrix 

(POUT = UIN IS permitted) . 


N 


Dimension of the matrices involved. 
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20. UD2SIG (tl— D factors to sxgmas) 

Purpose 

To compute variances from the U-D factors of a matrix. 

f~~ 

^ CALL UD2SIG(U,N,SIG,TEXT,NCT) 


Argument Definitions 
U(N*(H+l)/2) 


N 

SIG(N) 
TEXT ( ) 


NCT 


Input vector stored array containing 
the U-D factors. The D (diagonal) 
elements are stored on the diagonal 
of U. 

Dimension of the U matrix 

Output vector of standard deviations 

Output label of field data characters, 
which precedes the printed vector of 
standard deviations . 

Number of characters of text, 
O.LE.NCT.LE.126. If NCT = 0, no 
sigmas are printed, i.e. nothing is 
printed. 


Functional Description 

If U and D are written as doubly subscripted matrices then 

SIG(J) = (D(J,J) + D(K,K) [U(J,K)]^ ) 

^ K=J+1 / 


If NCT.GT.O a title is printed, followed by the sigmas. 
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21, UTINV (Upper triangular matrix inverse) 

Purpose 

To invert an upper triangular vector stored matrix and store the 
result in vector form. The algorithm is so arranged that the result 
can overwrite the input. 

CALL UTIN V ( RIN , N , ROUT ) 


Argument Definitions 

RIN(l'i*(S+l) /2) Input vector stored upper triangular 

matrix 

N Matrix dimension 

R0UT(N*(N+1) /2) ^ Output vector stored upper triangular 

matrix inverse (ROUT = RIN is per- 
mitted 

Remarks and Restrictions 

111 conditioning is not tested, but for nonsingular systems the 

result is as accurate as is the full rank singiilar value decomposition 

inverse. Singularity occurs if a diagonal is zero. The subroutine 
terminates when it reaches a zero diagonal. The columns to the left 
of the zero diagonal are, however, inverted and the result stored 
in ROUT. 

This routine can also be used to produce the solution to RX = Z. 
Place Z in column N+1 (viz, RIN (N* (N+1) /2+1) = Z(l), etc.), define 
RIN((N+1)XN+2)/2) - -1 and call the subroutine using N+1 instead of 
N. On return the first N entries of column N+1 contain the solution 
(e.g. ROUT (N* (N+1) /2+1) = X(l) , etc.). 

Because matrix inversion is numerically sensitive we recommend 
using this subroutine only in double precision. 
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Functional Description, 

The matrix inversion is accomplished using the standard back 
substitution method for inverting triangular matrices, cf. the book 
references by Lawson and Hanson, [1] or .•iierman [3]. 
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22. UTIROW (Upper triangular inverse, inverting only the upper rows) 
Purpose 

To compute the inverse of a vector stored upper triangular 
matrix, when the lower right comer triangular inverse is given. 

CALL UTIROW(RIN,N,ROUT,HRY) 


Argument Definitions 
RIN(N*(N+1) /2) 


N 

R0UT(N*(N+1) /2) 


NRY 


Remarks and Restrictions 


Input vector stored upper triangular 
matrix. Only the first N ~ NRY rows 
are altered by the algorithm. 

I 

Matrix dimension. 

Output vector stored upper triangular 
matrix inverse. On input the lower 
NRY dimensional right comer contains 
the given (known) inverse. This lower 
right corner matrix is left unchanged. 
(ROUT = RIN is permitted.) 

Number of rows, starting at the bottom, 
that are assumed already inverted. 


The purpose of this subroutine is to complete the coriputation 

of an upper triangular matrix inverse, given that the lower right 

corner has already been inverted. Part of the input, the rows to 

be inverted,' are inserted via th.e matrix RIN. The portion of the 

matrix that has already been inverted is entered via the jjiatrix ROUT. 
It may seem odd that part of the Input matrix is put into RIN and 

part into ROUT. The reasoning behind this decision is that RIN 

represents the input matrix to he inverted (it just happens that 

we do not make use of the lower right triangular entries) ^ ROUT 

represents the inversion result, and therefore that portion of the 

inversion that is given should be entered in this array. 
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111 conditioning is not tested, but for nonsingular systems the 

result is accurate. Singularity halts the algorithm if any of the 

first N-NRY diagonal elements is zero. If the first zero encountered 

moving up the diagonal (starting at N-NRY) is at diagonal j then the 

rows below this element will be correctly represented in ROUT. 

To generate estimates do the following: put N+1 into the matrix 

dimension argument; in the first N-NRY rows of the last column of 

RIN put the right hand side elements of the equation R^x +.R^^y = z^ 

(i.e., R , R , and z make up the first N-NRY rows of RIN); in the 

next NRY entries of ROUT, beginning in the (N-NRY+l)st element, put 

y (i.e., R ^ and y make up rows N— NRY+1, . . , ,N of ROUT); and 
est y est 

ROUT ( (N+1) (N+2)/2) = -1. On output, the last column of ROUT will 


contain x y ^ and -1. 
est est 

When NRY = 0 this algorithm is equivalent to subroutine UTIiW. 


Functional Description 

The matrix inversion is accomplished using the standard back 
substitution method. The computations are arranged row-wise, starting 
at the bottom (from row N-NRY, since it is asstamed that the last NRY 
rows have already been inverted). 
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23. WGS (Weighted Gram-Schmidt matrix triangularization) 

Purpose 

To compute a vector stored U-D array from an input rectangular 

T T 

matrix W, and a diagonal matrix D so that W D W = UDU . 

W XV 

GALL WGS(W,IMAXI^,IW,JW,DW,U,V) 


Argument Definitions 
W(IW,J¥) 

IMAXW 

DW(JW) 

U(IW*(IW+l)/2) 

V(JW) 


Input rectangular matrix, destroyed by 
the computations 

Roxv dimension of input W -matrix, 

IMAXW. GS.IW 

Diagonal input matrix; the entries 
are assumed to be nonnegative. This 
vector is unaltered by the computations 

Vector stored output D-D array 

Work vector in the computation 


Remarks and Restrictions 

The algorithm is not numerically stable when negative DW weights 

are used; negative weights ar^, however, allowed. If JW is less than 
IW (more rows than columns), the output U-D array is singular; with 
IW-JW zero diagonal entries in the output U array. 

Functional Description 

A D^-orthogonal set of row vectors, <f>^, ‘^’xw’ con- 

T 

strutted from the input rows of the W matrix, i.e., W=U<t), , ^D^<ji = D. 

The construction is accomplished using the modified Gram-Schiaxdt 
orthogonal construction (see refs. [1] or [3]). This algorithm is 

reputed to have excellent numerical properties. Note that the (f> 


vectors are not of interest in this routine, and they are overwritten; 

The V vector used in the program houses vector IW-j+l of <j) at step 3 of 

algorithm. The fact that the computed <}> vectors may not be D orthogonal 
is of no import in regard to the U and D computed results. 


62 



77-26 


V. FORTRAN Subroutine Listings 
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C 


C 


C 


c 


subroutine A6TRN (UINrUOUTrNr C r V) 

A6EE-TURNER U“D FACTOR RANK 1 UPDATE 


( UOUT ) +DOUT+ ( UOUT ) **T= (UIN ) *DIN* ( UIN ) ♦♦T+C*V+V**T 


UlN(N*(N+l)/2) 

U0UT{N+(N+l)/2) 

N 

C 

V(N) 


INPUT' vector STORED POSITIVE SEMI-DEFINITE U-D 

array» with q elements stored on the diagonal 

OUTPUT VECTOR STORED POSITIVE (POSSIBLY) SEmI- 
DEFINITE U-D result. UoUT=UIN is PERMITTED 

dimension OF the state 

SCALAR, should BE NON-nEGATIVE 

C IS destroyed during the process 

input vector for rank one modification. V IS 

destroyed during the process 


cognizant persons: G.J.BIERMAN/M.W.NEAD (JPL» FEB, 1977) 


IMPLICIT DOUBLE PRECISION (A-H»0-2) 
dimension UIN(1)» U0UT(1)» V(1) 


Z=0.0 

IF (C.EQ.Z) return 

JJ=N*(N+1)/2 
DO 50 J=N»2»“1 
S=V(J) 

D=UIN(JJ)+C*S*S 
IF (D) 5rl0f30 
5 WRITE (6»100) 

return 

10 

WRITE (6»110) 

DO 20 K=1»J 
20 UOUT(JJ+K)=Z 
GO TO 50 
30 B=C/D 

8ETA=S*B 

C=B*UIN(JJ) 

UOUT( JJ)=D 

JM1=J-1 
DO 40 I=1»JM1 

V ( I ) =V ( I ) -S*UIN ( JJ+I ) 

40 UOUT { JJ+I ) =UIN ( JJ+I ) +BETA*V ( I ) 

50 CONTINUE 

U0UT(1)=UIN(1)+C*V(1)**2 

RETURN 

100 format (1H0»10X»** * * ERROR RETURN DUE TO A COMPUTED NEGATIVE 
IPUTED DIAGONAL IN AGTRN * ♦ *f) 

110 format (lH0fl0X»** * ♦ note: U-D RESULT IS SINGULAR ♦ * ♦») 

end 


agtrnoio 

AGTRN020 

AGTRN030 

AGTRN040 

AGTRN050 

AGTRN060 

AGTRN070 

AGTRN080 

AGTRN090 

AGTRNIOO 

AGTRNIIO 

A6TRN120 

AGTRN130 

AGTRN140 

AGTRN150 

AGTRN160 

AGTRN170 

AGTRN180 

AGTRN190 

AGTRN200 

AGTRN210 

AGTRN220 

AGTRN230 

AGTRN240 

AGTRN250 

A6TRN260 

AGTRN270 

AGTRN280 

AGTRN290 

AGTRN300 

AGTRN310 

AGTRN320 

AGTRN330 

AGTRN340 

AGTRN350 

AGTRN360 

A6TRN370 

AGTRN380 

AGTRN390 

AGTRN400 

A6TRN410 

AGTRN420 

AGTRN430 

AGTRN440 

AGTRN450 

AGTRN460 

AGTRN470 

AGTRN480 

A6TRN490 

AGTRN500 

COMAGTRN510 

AGTRN520 

AGTRN530 

AGTRN540 
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60 

70 

80 

90 

100 


SUBROUTINE A2A1 ( A» lA » IR»LA » NAMA» A1 » lAl »LA1 » NAMAl) 

SUBROUTINE TO REARRANGE THE COLUMNS OF A(lR.LA)f IN NAMA ORDER 
AND PUT THE RESULT IN AKlR.LAll IN NAMAI ORDER. ZERO COLUMNS 
ARE inserted IN Al CORRESPONDING TO THE NEWLY DEFINED NAMES, 

A(iR»LA) Input rectangular matrix 

lA Row DIMENSION OF A» IR.LE.lA 

IR NO. OF ROWS Of A THAT ARE TO BE REARRANGED 

LA NO. OF PARAMETER NAMES ASSOCIATED WITH A 

NAMA (LA) PARAMETER NAMES ASSOCIATED WITH A 

AKIRfLAl) OUTPUT RECTANGULAR MATRIX 

A AND Al cannot SHARE COMMON STORAGE 
lAl ROW DIMENSION OF Air IR.LE.lAl 

LaI no. of parameter names associated WITH Al 

NAMAI (LAI) INPUT LIST OF PARAMETER NAMES TO BE ASSOCIATED 
WITH THE OUTPUT MATRIX Al 

COGNIZANT persons: G. J .BIERMAN/M. W.NEAD (JPL* sept. 1976) 

DIMENSION ACIArDr NAMA(i)r A 1 { I Al r 1 ) r NAMAI ( 1 ) 
implicit DOUBLE PRECISION (A-HrO-Z) 

ZERO=0. 

DO 100 J=lrLAl 
DO 60 I=ifLA 

IF (NAMA(I) .EQ.NAMAl(J) ) GO TO 60 
CONTINUE 


DO 70 K=lfIR 
Al (Kr J)=ZERO 
60 TO 100 
00 90 K=lrlR 
AlCKr J)=A(K>D 
CONTINUE 

RETURN 

END 


0 ZERO COL, corrfs. To new name 
(3 COPY COL. ASSOC. WITH OLD NAME 


A2A10010 

A2A10020 

A2A10030 

A2A10040 

A2A10050 

A2A10060 

A2A10070 

A2A10080 

A2A10090 

A2A10100 

A2A10110 

A2A10120 

A2A10130 

A2A10140 

A2A10150 

A2A10160 

A2A10170 

A2A10180 

A2A10190 

A2A10200 

A2A10210 

A2A10220 

A2A10230 

A2A10240 

A2A10250 

A2A10260 

A2A10270 

A2A10280 

A2A10290 

A2A10300 

A2A10310 

A2A10320 

A2A10330 

A2A10340 

A2A10350 

A2A10360 
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SUBROUTINE COMBO (R»L1 » NAMI »L? . NAM2f Ar lA ► LA.NAMA ) 


TO rearrange a vector stored triangular matrix and store 

THE RESULT IN MATRIX A. THE DIFFERENCE BETWEEN THIS SUB- 
ROUTINE AND R2A IS that THFRE THE NAMELIST FOR A IS INPUT. 
HERE IT IS determined BY COMBINING THE LIST FOR R WITH 
A LIST OF desired NAMES. 


R(Ll*(Ll+l)/2) 

LI 

NAMKLl) 

L2 

NAM2(L2> 


A(L1»LA) 


lA 

LA 


NAMA(LA) 


INPUT VECTOR STORED UPPER TRIANGULAR MATRIX 
NO. OF parameters IN R <AND IN NAMl) 

NAMES ASSodlATED WITH R 
NO. OF PARAMETERS IN NAM? 

PARAMETER NAMES THAT ARE TO BE COMBINED WITH R 
(NAMl LIST). THESE NAMES MAY OR MAY NOT RE IN 
NAMl. 

OUTPUT ARRAY CONTAINING THE REARRANGED 
R MATRiXf Ll.LE.IA. 

ROW DIMENSION OF A 

NO. OF PARAMETER NAMES In NAMA» AND THE 
COLUMN DIMENSION OF A. LA=L1+L2-N0. NAMES 
COMMON TO naMI AND NAM2. LA IS COMPUTED AND 
OUTPUT. 

parameter names associated with THE OUTPUT A 
MATRIX. consists OF NAMES IN NAMl NOT IN 
NAM2 followed BY NAM2, 


COGNIZANT persons: G. J.BIERMAN/M.W.NEAD tJPLt SEPT. 197ft) 

IMPLICIT DOUBLE PRECISION (A-HrO-2) 

dimension R(l)r A(IA»1)» NAm 1(1)» NAM2{l)r NAMA(l) 

C 

zero=o.o 

K=1 

DO 100 I=1»L1 
DO 50 J=1»L2 

IF (NAMl(I) .EQ.NAM2(J) ) GO TO 100 
50 CONTINUE 

NAMA(K)=NAM1(I) 

JJ=I+(I-l)/2 
DO 60 L=lfl 
60 A(L»K)=R(JJ+L) 

IF (I.EQ.Ll) GO TO flO 
IPI - I+l 
DO 70 L=IP1»L1 
70 A(L.K) = ZERO 
80 k=K+l 

100 continue 

C NAMES unique TO NAMl ARE NOW IN NAMA 

do 200 J=1»L2 
DO ISO I=lfLl 

IF {NAM2(J) .EQ.NAMKD ) GO TO 170 
150 CONTINUE 

NAMA(K)=NAM2( J) 

DO 160 L=1»L1 
160 A(L»K)=ZERO 


COMBOOlO 

COMBO020 

COMB0030 

COMR0040 

COMB0050 

COMB0060 

COMB0070 

COMB0080 

COMB0090 

COMBOlOO 

COMBOllO 

COMB0120 

COMB0130 

COMB0140 

COMB0150 

COMB0160 

C0M8017Q 

COMBO180 

COMB0190 

COMBO20Q 

COMBO210 

COMB0220 

COMB0230 

COMB0240 

COMBO250 

COMR0260 

COMB0270 

C0MB0280 

COMR0290 

COMBO300 

COMB0310 

COMB0320 

COMBO330 

COM8O340 

C0MB035Q 

COMB0360 

COMB0370 

COMBO380 

COMB0390 

COMR0400 

COMBO«J10 

COMB0420 

C0MBO430 

COMB0440 

COMB0450 

COM80460 

COMB0470 

COMBO480 

COMB0490 

COMB0500 

COMR0510 

COMB0520 

COMB0530 
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C NAMES unique TO NAM? ARE 

60 TO 190 

170 NAMA<K)=NAM2C J) 

C ■ ‘ LOCATE diagonal OF PRECEDIN6 

JJ=I4:(I-1)/2 

DO 180 L=1»I ' 

180 A(L»K)=R'CJJ+L) 

IF (I.EQ.LI) go TO 190 
IP1=I+1 

DO 185 L=IP1»L1 
185 A(LrK)=ZERO" 

190 ' K=K+1 

200 CONTINUE 
llA=K-l 

C ‘ ' NAMES MUTUAL TO NAM]' AND 

RETURN ' 

END 


NOW IN NAMA COMB0540 

COMB0550 

COMB0560 

column COMB0570 

COMB0580 

COMR0590 

COMB0600 

COMR0610 

COM80620 

COMBO630 

COM80640 

COMR0650 

COMB0660 

COMB0670 

NAM2 ARE NOW IN NAMA COMB06B0 

COMR0690 

COMR0700 
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SUBROUTINE C0V2Rl(UfN) 

TO CONSTRUCT THE UPPER TRIANGULAR ChOLESKY FACTOR OF A 
POSITIVE SEMI-DEFINITE MATRIX. BOTH THE INPUT COVARIANCE 
AND THE OUTPUT CHOLESKY FACTOR (SQUARE ROOT) ARE VECTOR 
STORED. THE OUTPUT OVERWRITES THE INPUT, 

COVARIANCE { INPUT) =U+U**T (U IS OUTPUT). 

IF THE INPUT covariance IS SINGULAR THE OUTPUT FACTOR HAS 
ZERO COLUMNS. 

U(N*(N+l)/2) contains THE INPUT VECTOR STORED COVARIANCE 

MATRIX (ASSUMED POSITIVE DEFINITE) AND ON OUTPUT 
IT contains the upper TRIANGULAR SQUARE ROOT 
FACTOR. 

N DIMENSION OF THE MATRICES INVOLVED 

COGNIZANT persons: G. J.BIFRMAN/M. W.NEAD (JPL» FEB. 1977) 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) 

DIMENSION U(l) 

Z£RO=0.0 

ONE=l. 

JJ=N^=(N+l)/2 

JJN=JJ 

DO 5 J=N»2»-I 

IF (U(JJ) .LT.ZERO) U(JJ)=ZFR0 
U(JJ)= SQRT(U(JJ)) 

IF (U(JJ) .GT.ZERO) ALPHA=0NE/U(JJ) 


KK=0 

JJN=JJ-J 

JM1=J-1 

DO 4 K=l.JMl 

U ( J JN+K ) =ALPHA+U { J JN+K ) 
S=U(JJN+K) 

DO 3 I=1»K 

3 U(KK+I)=U(KK+I)-S*U(JJN+I) 

4 KK=KK+K 

5 JJ=JJN 

IF (U(l) .LT.ZERO) U(1)=ZER0 
UU)= SQRT(Ud)) 

REtURN 

END 


Q NEXT diagonal 

0 JJN+K=(K»sJ) 
ra KK+I=(I»K) 


COV2R010 

COV2R020 

COV2R030 

COV2R040 

COV2R050 

COV2R060 

COV2R070 

COV2R080 

COV2R090 

COV2R100 

COV2R110 

COV2R120 

COV2R130 

COV2R140 

COV2R150 

COV2R160 

COV2R170 

COV2R1BO 

COV2R190 

COV2R200 

COV2R210 

COV2R220 

COV2R230 

COV2R240 

COV2R2SO 

C0V2R260 

COV2R270 

COV2R280 

COV2R290 

C0V2R300 

COV2R310 

COV2R320 

COV2R330 

COV2R340 

COV2R350 

C0V2R360 

COV2R370 

COV2R30O 

COV2R390 

COV2R400 

COV2R410 

COV2R420 

COV2R430 

COV2R440 

COV2R450 

COV2R4GO 
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SUBROUTINE C0V2UD (U»N) 

COV2U010 
COV2U020 
COV2U030 
COV2U040 
COV2U050 
COV2U060 
COV2U070 
COV2U080 
COV2U090 
COV2U100 

SINGULAR INPUT COVARIANCES RESULT IN OUTPUT MATRICES WITH ZERO COV2UUO 



COLUMNS 

COV2U120 

COV2U130 

COV2U140 


COGNIZANT persons: G.J.BIFRMAN/R.A.JACOrSON (JPL» FEB. 1977) 

C0V2U150 

COV2U160 


IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

COV2U170 

COV2U180 


DIMENSION U(l) 

COV2U190 

COV2U200 


z=o.o 

COV2U210 


oNe=i.o 

COV2U220 

COV2U230 


jJ=N*(N+l)/2 

CnV2U240 


DO 50 vJ=Nr2r-l 

COV2U250 


ALPHA=Z 

COV2U260 


IF (U(JJ).LT.Z) U(JJ)=Z 

COV2U270 


IF (U(JJ).GT.Z) ALPHA=ONE/U( JJ) 

COV2U280 


JU— vJU— u 

COV2U290 


KK=0 

COV2U300 


KJ=JJ 

COV2U310 



COV2U320 


DO 40 K=1»JM1 

COV2U330 


KU=KJ+1 

COV2U340 


BETA=U(KJ) 

COV2U350 


U(KJ)=ALPHA*U(KJ) 

COV2U360 


IJ=JJ 

COV2U370 


IK=KK 

COV2U380 


DO 30 1=1. K 

COV2U390 


IK=IK+1 

COV2U400 


IJ=IJ+1 

COV2U410 

30 

U(IK)=U(IK)-BETA*U(IJ) 

COV2U420 

40 

KK=KK-t-K 

COV2U430 

50 

CONTINUE 

COV2U440 


IF {U(l).LT.Z) U(l)=7 

COV2U450 


RETURN 

COV2U460 


end 

COV2U470 


TO OBTAIN THE U-D FACTORS OF A POSITIVE SEMI-DEFTNITE MATRIX. 
THE INPUT matrix VECTOR STORED IS OVERWRITTEN BY THE OUTPUT 
U-D Factors which are also vector stored. 

U{N*(N+1)/2) contains input vector stored covariance matrix. 

ON OUTPUT IT contains ThE VECTOR STORED U-D 
COVARIANCE FACTORS. 

N MATRIX DIMENSION 
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SUBROUTINE C2C (C»IC»Ll»NAMi.L?rNAM2) 

SUBROUTINE TO REARRANGE THE ROWS AND COLUMNS OF MATRIX 
C(Li»Ll) IN NaMI order 'and PUT THE RESULT IN 
C(L2»L2) IN NAM2 ORDER. ZERO COLUMNS AND ROWS ARE 
ASSOCIATED WITH OUTPUT DEFINED NAMES THAT ARE NOT CONTAINED 
IN NAMI. 

C(LlfLl) INPUT MATRIX 

IC ROW DIMENSION OF Ct IC .GE.L=MAX (LI »L2) 

LI NO. OF PARAMETER NAMES ASSOCIATED WITH THE INPUT ( 

NAMl(L) PARAMETER NAMES ASSOCIATED WITH C ON INPUT, (ONLY 

The FIRST LI entries APPLY TO THE INPUT C) 

L2 NO. OF parameter NAMES ASSOCIATED WITH THE OUTPUT 

NAM2(L2) PARAMETER NAMES ASSOCIATED WITH THE OUTPUT C 

COGNIZANT PERSONS! G, J.BIERMAN/M, W.NEAD (JPL» SEPT, 1976) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION C(ICrl)» NAMl(l)» NAM2(1> 

Z£RO=0, 

L=MAX(LlrL2) 

IF (L.LE.Ll) GO TO 5 
NM=L1+1 
DO 1 K=NM*L 
1 NAMl(K)= ZFRO 
5 DO 90 J=1»L2 
DO 10 I=1»L 

IF (NAMl(I) .EQ.NAM2(J) ) 

10 CONTINUE 

60 TO 90 

30 IF (I.EQ.J) GO TO 90 
DO 40 K=lrL 
H=C(Kf J) 

C(Kf J)=C(K»I) 

40 C(KrI)=H 

DO 80 K=lrL 
H=C(JfK) 

C(J»K)=C(I,K) 

80 C(I»K)=H 

NM=NAM1(I) 

NAMl(I)=NAMl(J) 

NAMI ( J)=NM 
90 CONTINUE 

FIND NAM2 NAMES NOT IN NAMI AND SET CORRESPONDING ROWS AND 
COLUMNS To ZERO 

DO 120 J=1»L2 
DO 100 I=1»L 

IF (NAMl(I) .EO.NAM2(J) ) GO TO 120 
100 CONTINUE 

DO 110 K=1,L2 
C(J»K)=ZFRO 
110 C(KrJ)=ZERO 

120 CONTINUE 


ra ZERO REMAINING NAMI LOCNS 


GO TO 30 


Q interchange columns I AND J 


a interchange rows i and j 


0 INTERCHANGE LABELS I AND J 


return 

end 
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C2C00010 
C2CO,0020 
C2C00030 
C2C00040 
C2C00050 
CPC00060 
C2C00070 
C?C00080 
C2C00090 
: C2C00100 
C2C00110 
C2C00120 
CC2C00130 
C2C00140 
C2C00150 
C2C00160 
C2C00170 
C2C00180 
C9C00190 
C2C00200 
C2C00210 
C2C00220 
C2C00230 
C2C00240 
C2C00250 
C2C00260 
C2C00270 
C2C00280 
C2C00290 
C2C00300 
C2C00310 
C2C00320 
C2C00330 
C2C00340 
C2C00350 
C2C00360 
C2C00370 
C2C00380 
C2C00390 
C2C00400 
C2C00410 
C2C00420 
C2C00430 
C2C00440 
C2C00450 
C2C00460 
C2C00470 
C2CQ0480 
C2C00490 
C2C00500 
C2C0051O 
C2C00520 
C2C00530 
C2C00540 
C2C00550 
C2C00560 
C2C00570 
C2C00580 
C2C00590 
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C 

C 


SUBROUTINE IMF2R (PfN) 

TO CHOLESKY FACTOR AN INFORMATION MATRIX 

COMPUTES A LOWER TRIANGULAR VECTOR STORED CHOLESKY FACTORIZATION 
OF A POSITIVE SEMI-DFFINITE MATRIX. P=R(**T)R» R UPPER TRIANGULAR 
BOTH MATRICES ARE VECTOR STORED AND THE RESULTS OVERWRITES 
THE INPUT 


PtN*(N+l)/2) 


N 


10 


ON INPUT THIS IS A POSITIVE SfMI-OEFINITE MATRiXf 
AND ON OUTPUT IT IS A TRIANGULAR FACTOR. IF THE 
input matrix Is singular The output matrix will 
HAVE ZERO diagonal ENTRIES 
DIMENSION OF MATRICES INVOLVED 


cognizant person: G.J.BIERKAN/M.W.NEAD 
implicit double precision (A-HrO-Z) 
DIMENSION P<1) 


Z=0,0 

0NE=1.0 

JJ=0 

NN=N*tN+l)/2 

NM1=N~1 

DO 10 J=lrNMl 

IF (P(JJ).LT.Z) P(JJ3=Z 
P(JJ)=SQRT(P(JJ) ) 

ALPHA=Z 

IF (P(JJ).GT.Z) ALPHA=ONE/P(JJ) 

JK=NN+J 

JP1=J+1 

J1S=JK 

DO 10 K=NrJPl»~l 
JK=JK-K 

p{JK)=ALPHA*PUK) 

8ETA=P{JK) 

KI=NN+K 

JI=JIS 

DO 10 I=N»K»-1 
KI=KI~I 
JI=JI-I 

P(KI)=P(KI)~P( JD+BETA 


(JPL»FEB.1977) 


Q UJ— ( J r U ) 

0 JK=(J.K) 

Q JI.S=(J»I) START 


IF (P(NN).LT.Z) PtNN)=Z 
P<NN3=SQRT(P(NN) ) 

RETURN 

END 


INF2R010 

■INF2R020 

INF2R030 

INF2R040 

INF2R050 

.INF2R060 

INF2R070 

INF2R080 

INF2R090 

INF2R100 

INF2R110 

INF2R120 

INF2R130 

INF2R140 

INF2R150 

INF2R160 

INF2R170 

INF2R180 

INF2R190 

INF2R200 

INF2R210 

INF2R220 

INF2R230 

INF2R240 

INF2R250 

1NF2R260 

INF2R270 

INF2R2B0 

IMF2R290 

INF2R300 

INF2R310 

INF2R320 

IMF2R330 

IMF2R340 

INF2R350 

INF2R360 

INF2R370 

IMF2R380 

IMF2R390 

INF2R400 

IMF2R410 

INF2R420 

INF2R430 

INF2R440 

IMF2R450 

INF2R460 

INF2R470 

INF2R48C 

IMF2R490 

INF2R50C 
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of ™ 


SUBI<OUTINE PERMUT (A»IA» IR»L1 ►NAMa»L2rNAM23 

SUBROUTINE TO REARRANGE PARAMETERS OF A(TRrLl>i NAMl ORDER 
TO A(lRrL2)f NAM2 ORDER. ZERO COLUMNS ARE INSERTED 
corresponding to the newly defined NAMES, 


AURfL) 

lA 

IR 

Li 

NAMl(L) 

L2 

NAM2 


INPUT RECTANGULAR MATRIX* L=MAX(Li»L2) 

ROW DIMENSION OF A* lA.GE.lR 

NUMBER OF ROWS OF A THAT ARE TO BE REARRANGED 
NUMBER OF PARAMETER NAMES ASSOCIATED WITH THE INPUT 
A matrix 

PARAMETER NAMES ASSOCIATED WITH A ON INPUT 
(ONLY THE FIRST LI ENTRIES APPLY TO THE INPUT A) 
number OF parameter NAMES ASSOCIATED WITH THE OUTPUT 
A matrix 

parameter 


names associated with the OUTPUT A 


COGNIZANT persons; G. J.BIERMAN/M. W.NEAD (JPLr SEPT, 1976) 

IMPLICIT DOUBLE PRECISION (A-H*0-Z) 

DIMENSION A(IA»1)» NAMKD* NAM2(1) 

2ERO=0, 

L=MAX(L1»L2) 

IF (L.LE.Ll) GO TO 50 
NM=L1+1 
DO 40 K=NM»L 
40 NAM1(K)=0 

50 DO 100 J=lfL2 
DO 60 I=1*L 

IF (NAMid) .EQ.NAM2(J) ) 

60 CONTINUE 

GO TO 100 
65 CONTINUE 

IF (I.EQ.J) GO TO 
DO 70 K=1*1R 
W=A(K*J) 

A(K» J)=A(K.I) 

70 A(Krl)=W 

NM=NAM1(1) 

NAM1(I)=NAM1(J) 

NAMl ( J)=NM 
100 CONTINUE 

REPEAT TO FILL NEW COLS 

DO 200 J=1»L2 
DO 160 1=1, L 

IF (NAMl(I) ,E0,NAM2(J)) gO TO 200 
160 CONTINUE 

DO 170 K=lrlR 
170 A(K*U)=ZERO 

200 CONTINUE 

return 
end 


0 ZERO remaining NAMl LOCS 


GO TO 65 


100 

Q interchange cols I AND J 


Q interchange I and J COL. LABELS 


PERMUOlO 

PERMU020 

PERMU030 

PERMU040 

PERMU050 

PERMU060 

PFRMU070 

PERMU080 

PERMU090 

PERMUlOO 

PERMUilO 

PERMU120 

PERMU130 

PFRMU140 

PERMU150 

PERMU160 

PERMU170 

PFRMU180 

PERMU190 

PERMU200 

PERMU210 

PFRMU220 

PFRMU230 

PERMU240 

PFRMU250 

PERMU260 

PERMU270 

PFRMU280 

PERMU290 

PERMU300 

PFRMU310 

PERMU320 

PERMU330 

PERMU340 

PERMU350 

PERMU360 

PERMU370 

PERMU380 

PERMU390 

PFRMU400 

PERMU410 

PFRMU420 

PFRMU430 

PFRMU440 

PERMU450 

PERMU460 

PERMU470 

PERMU480 

PERMU490 

PERMU500 

PFRMU510 

PFRMU520 

PERMU530 
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C 


c 


c 


c 


c 

c 


c 


SUBROUTINE RlNCON XRIN»N»R0UT*CNB) 

TO COMPUTE THE INVERSE OF THE UPPER TRIANGULAR VECTOR STORED 
INPUT MATRIX RIN AND STORE The RESULT IN ROUT. (RlNrROUT IS 
PERMITTED) AND TO COMPUTE A CONDITION NUMBER ESTIMATE. 
CNB=FR08 .norm ( R ) *FROR . NORM ( R**- l ) . 

THE FROBENIUS NORM IS THE SQUARE ROOT oF THF SUM OF SQUARES 
OF THE ELEMENTS. THIS CONDITION NUMBER BOUND IS USED AS 
an UPPER BOUND AND IT ACTS AS A LOWER ROUND' ON THE ACTUAL 
CONDITION number OF THE PROBLEM. (SEE THE BOOK ’SOLVING LEAST 
SQUARES’* BY LAWSON AND HANSON) 

RlN(N*(N+l)/2) INPUT VECTOR STORED UPPER TRIANGULAR MATRIX 
N dimension of r MATRICES 

R0UT(N=«'(N+l)/2) OUTPUT VECTOR STORED UPPER TRIANGULAR MATRIX 
inverse (RIN=RoUT is PERMITTED) 

cnb condition numbfr bound 

cognizant persons: G.J.BIERMAN/M.W.NEAD (JPLrFEB. 1977) 

SUBROUTINES REQUIRED: UTINV 

IMPLICIT DOUBLE PRECISION (A-H*0-Z) 

DIMENSION RIN(D* ROUT(l) 


2 = 0.0 

NT0T=N+(N+1)/2 

RNM=Z 

DO 10 J=1*NT0T 
10 RNM=RNM+RIN(J)**2 

CALL UTINV (RlNrN*ROUT) 

RNMOUT=Z 

DO 20 J=1*NT0T 

20 RNM0UT=RNM0UT+R0UT< J)**2 

CNB=SQRT (RNM*RNM0UT) 

WRITE (6*30) CN8 
RETURN 


RTNCOOlO 

RINC0020 

RINC0030 

RTNC0040 

RINC0050 

R1NC0060 

RINC0070 

RTNC0080 

RINC0090 

RiNCOlOO 

RINCOIIO 

RINCO120 

R1NCO130 

RINCO140 

RINCO150 

R1NCO160 

RINCO170 

RTNC0180 

RINC0190 

RTNC0200 

RTNCO210 

RTNCO220 

RINCO230 

RTNCO240 

RINC0250 

RINC0260 

PINC0270 

RINCO280 

RINC0290 

RINC0300 

R1NCO310 

RINCO320 

RINC0330 

RINCO340 

RINCO350 

RINCO360 

RINC0370 

RINC0380 

RTNC0390 

RINC0400 

RINCO410 

RINCO420 


30 FORMAT(1HO*5X* ’CONDITION NUMBFR 
IION NUMBER. LE. cnb’ */) 

END 


BOUND=' *D18.10*2X* *CNB/N.LE.CONniTRlNC0430 

RINC0440 

RTNC0450 
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C 

C 


c 

c 


1 

2 


5 

10 


SUBROUTINE RI2C0V (RTNVr N»SIS»C0V0UT»K0V) 

RI2C0010 

TO COMPUTE the COVARIANCE MATRIX AND/OR THE STANDARD DEVTATIONSRI2C0020 
OF A VECTOR STORED UPPER TRIANGULAR SQUARE ROOT COVARIANCE RI2C0030 
MATRIX. THE OUTPUT COVARIANCE MATRIX IS ALSO VECTOR STORED. RI2C0040 

RI2C0050 

RlNV(N*(N+l)/2) INPUT VECTOR STORED UPPER TRIANGULAR COVARI- R12C0060 

ANCE square ROOT. (RINV=R INVERSE IS THE RI2C0070 
inverse of the SRIF MATRIX) RI2C0080 

N dimension OF THE RINV MATRIX RI2C0090 

SIG(N) OUTPUT vector OF STANDARD DEVIATIONS RT2C0100 

C0V0UT(N*(N+l)/2) OUTPUT VFCTOR STORED COVARIANCE MATRIX RI2C0110 

(COVOUT = RINV IS ALLOWED) RI2C0120 

KOV .GT.O COMPUTE COVARIANCE AND SIGMAS USING KOv ROWS RI2CO130 


OF RiNV. 


RT2C0140 

.LT.O COMPUTE ONLY THE SIGMAS USING 

KOV ROWS OF 

RI2CO150 

RINV. 


RI2C0160 

.EO.O NO covariance* BUT ALL SIGMAS 

(E.G. USE 

RI2C0170 

N ROWS OF RINV). 


RT2CO180 

RI2C0190 

COGNIZANT persons: G. J.SIERMAN/M. W . NEAD (JPL» sept, 1976) 

RI2C0200 

RT2C0210 

IMPLICIT DOUBLE PRECISION CA-H»0-Z) 


RI2C0220 

DIMENSION RINV(l)* SIG(i)* COVOUTU) 


RI2C0230 

RI2C0240 

Z£RO=0,0 


RI2C0250 

LIM=N 


RI2C0260 

IF (KOV.NE.O) LlM=IARS(KOV) 


RI2CO270 

♦♦♦ COMPUTE SIGMAS 


RT2C0280 

IKS=0 


RI2C0290 

DO 2 J=1*LIM 


RT2C0300 

IKS=IKS+J 


RI2C0310 

5UM=ZER0 


RI2C0320 

IK=IKS 


RI2C0330 

DO 1 K=J*N 


RI2CO340 

SUM=SUM+R INV ( IK ) **2 


PI2C0350 

IK=IK+K 


PT2CO360 

SIG(J)=SQRT(SUm) 


PI2C0370 

RI2C0380 

IF (KOV.LE.O) RETURN 


R12C0390 

*** COMPUTE COVARIANCE 


RI2C0400 

JJ=0 


RI2C0410 

nmi=lim-i 


RI2CO420 

DO 10 J=lfNMl 


RI2C0430 

JJ=JJ+J 


RI2C0440 

COVOUT (JJ)=SIG(J)**2 


PI2C0450 

IJS=JJ+J 


RI2C0460 

JP1=J+1 


PI2CO470 

DO 10 I=JP1»N 


RI2C0480 

IK=IJS 


RI2C0490 

IMJ=I-J 


PI2C0500 

5UM=ZER0 


RI2C0510 

DO 5 K=I*N 


PI2C0520 

IJK=IK+IMJ 


RI‘2C0530 

SUM=SUM+r I NV ( IK ) *R I N V { UK ) 


PT2C0540 

IK=IK+K 


RI2C0550 

COVOUT ( US )=SUM 


RI2CO560 

IJS=IJS+1 


RT2C0570 

IF (KOV.EQ.N) C0V0UT(JJ+N)=SIG{N)**2 


RI2C0580 

RI2C0590 

RETURN 74 


RI2C0600 

END 


RI2C0610 
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subroutine R2A (R » LR f NAMR » A » I A » lA » NAMA ) 


TO PLACE THE TRIANGULAR VECTOR STORED MATRIX R INTO THE 
matrix A AND. TO ARRANGE THE COLUMNS' TO MATCH THE DESIRED 
NAMA PARAMETER LIST, .NAMES IN THE NAMA LIST THAT DO NOT 
CORRESPOND TO ANY NAME IN naMR HAVE ZERO ENTRIES IN THE 
CORRESPONDING A COLUMN. 


R(LR*(LR+1)>2) 

LR 

NAMR(L) 

A(IR»LA) 

lA 

LA 

NAMA (LA) 


input UPPER TRIANGULAR VECTOR STORED ARRaY 
DIMENSION OP R/ 

parameter names ASSOCIATED with R. ONLY THE 
FIRST LR entries APPLY TO Rr L=MAX (LR fLA) . 

matrix to house the rearranged r matrix 

ROW DIMENSION OF Ar lA.GF.LR 

NO, OF parameter names ASSOCIATED WITH ThE, 

OUTPUT A matrix 

parameter names for THE OUTPUT A MATRIX 


COGNIZANT persons: G. J.BIERMAN/M.W.NEAD (JPLr SEPT. 197ft) 


IMPLICIT DOUBLE PRECISION (A-H.O-Z) 
dimension R(1) »NaMR( 1) fA(lA»n .NAMA(l) 

ZERO=0. 

DO 5 J=ltLA 
DO 5 K=1»LR 

5 A(KfU)=2ER0 ' 0 ZERO A(LR»LA) 

DO 40 J=1»LA 
DO 10 I=1»LR 

IF (NAMR(I) .EQ.MAHA(J) ) GO TO 20 
10 CONTINUE 

GO TO 40 
20 JJ=I*(I-l)/2 
DO 30 K=1»I 
30 A(Kf J)=R( JJ+K) 

40 CONTINUE 

return 

END 


R2A00010 

RPA00020 

R2A00030 

P2A00040 

R2A00050 

P2A00060 

R2A00070 

P2A00080 

R2A00090 

R2A00100 

R2A00110 

R2A00120 

R2A00130 

R2A00140 

R2A00150 

R2A00160 

R2A00170 

RPA0D18O 

R2A00190 

R2A00200 

RPA00210 

R2A00220 

R2A00230 

R2A00240 

R2A00250 

R2A00260 

R2A00270 

R2A00280 

R2AO029O 

P2A00300 

R2A00310 

R2A00320 

R2A00330 

R2A00340 

R2A00350 

R2A00360 

R2A00370 

R2A00380 
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SUBROUTINE R2RA (RFNRfNAM»RA»NRA»NAMA) 


C 

C 


c 


c 


R?RA0010 

TO COPY THE UPPER LEFT (LOWpR RIGHT) PORTION OF A VECTOR . R2RA0020 

STORED UPPER triangular MATRIX R INTO THE LOWER RIGHT R2RA0030 

(UPPER LEFT) PORTION ,OF A VECTOR STORED TRIANGULAR R2RA0040 

MATRIX RA. ‘ R2RA0050 

R2RA0060 

R(NR*(NR+l)/2) INPUT VECTOR STORED UPPER TRIANGULAR MATRIX R2RA0070 

NR dimension OF R R2RA0080 

NAM(NR) NAMES ASSOCIATED WITH R RPRA0090 

. RA(NRA*(NRA+l)/2) OUTPUT VECTOR STORED UPPER TRIANGULAR MATRIX RpRAOlOO 

NRA DIMENSION ASSOCIATED WITH RA P2RA0110 

NAMA(NRA) NAMES ASSOCIATED WITH RA R2RA0120 

R2RA0130 

IF NRA=0 ON INPUT.# THEN NAMA(l) SHOULD HAVE THE FIRST NAME OF THE R2RA0140 
OUTPUT NAMELIST and THE NUMBFR OF NAMES IN .NAMA IS COMPUTED. R2RA0150 


THE LOWER RIGHT BLOCK OF R WILL RE THE UPPER LEFT BLOCK OF RA, ' R2RA0160 


IF NRA=LAST NAME OF THE UPPER LEFT BLOCK THaT IS TO BE MOVED# 
THEN THE UPPER BLOCK IS TO BE MOVED TO THE LOWER RIGHT POSITION. 
when USED IN THIS MODE NRA=NR ON OUTPUT. 

THE NAMES OF THE RELOCATED BLOCK ARE ALSO MOVED. THE RESULT 
CAN COINCIDE WITH R AND NAMA WITH NAM. 

COGNIZANT persons: G. J.BIERMAN/M.W.NEAD (JPL# SEPT. 197fi) 


IMPLICIT DOUBLE PRECISION (a-H#0-Z) 

DIMENSION R(l)#RA(i)# NAM(l)» NAMA(i) 

logical IS 

I S=, FALSE. 

LOCNSNAMAU) 

IS=FALSE CORRESPONDS TO MOVING UPPER LFT, CORNER OF R TO 
LOWER RT, CORNER OF RA 
IF (NRA.EQiO) GO TO 1 
LOCNsNRA 
ISs.TRUE. 

ISSTRUE CORRESPONDS TO MOVING LOWER LFT. CORNER OF R TO 
UPPER RT, CORNER OF RA 
1 DO 3 I=1»NR 

IF (NAM(I).FQ.LOCN) GO TO 4 

3 CONTINUE 

WRITE (6#100) 

100 FORMAT ,(1HG#20X» *NAMAa) NOT IN NAMELIST OF R MATRIX’) 
RETURN 

4 K=1 

KM1=K-1 

IF (IS) GO TO 15 

IJS=K*(K+l)/2-l 

NRA=NR-K+1 

IJA=Q 

KOLA=0 


R2RA0170 

R2RA0180 

R2RA0190' 

R2RA0200 

R2RA0210 

R2RA0220 

R2RA0230 

R2RA0240 

R2RA0250 

R2RA0260 

R2RA0270 

R2RA0280 

R2RA0290 

R2RA0300 

R2RA0310 

R2RA0320 

R2RA0330 

R2RA034Q 

R2RA0350 

R2RA0360 

R2RA0370 

P2RA0380 

R2RA0390 

R2RA0400 

R2RA0410 

R2RA0420 

R2RA0430 

R2RA0440 

R2RA0450 

R2RA0460 

P2RA0470 

R2RA048Q 

R2RA0490 

R2RA0500 

P2RA0510 

R2PA0520 

R2RA0530 

R2RA0540 
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DO 10 KOL=K»NR 


R?RA0550 


KOLA=KOLA+l 

NAMA (KOL-KMl \=NAM (KOL) 

DO -5 1R=1»'koLA ' 
IJA=IJA+1 

5 RA(IJA)=R(IJS+IR) 

10 IJS=IJS+KOL, 

return 

15 IJ=K*(K+l)/2 

lJA=NR*(NR+l)/2 

L=NR“KM1 

KOL=K 

DO 25 KOLA=NRfLf-l 
IJS=IJA 

NAMA(KOLA)=NAM(KOL) 
DO 20 IR=K0LA»L »-l 
RA(IdS)=R(IJ) 
IJS=IJS-1 

20 IJ=IJ-1 

IJA=IJA-KOLA 
25 K0L=K0L-1 

NRA=NR 

return 

END 


P2RA0560 

R2RA0570 

R2RA0580 

R2RA0590 

R2RA0600 

R2RA0610 

RPRA0620 

R2RA0630 

RRRA0640 

R2RA0650 

R2RA0660 

P2RA0670 

R2RA0680 

R2RA0690 

R2RA0700 

R2RA0710 

RRRA0720 

R2RA0730 

R2RA0740 

R2RA0750 

P2RA0760 

R2PA0770 

R2RA0780 

R2RA0790 

R2RA080D 
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C 


C 


C 


c 


SUBROUTINE RUDR (RIN »N» ROUT ► IS) 

FOR N.GT.O this subroutine TRANSFORMS AN UPPER TRIANGULAR VfCTOR 
stored SRIF matrix to U-D form# and when N.lT.O THE U-D VECTOR 
STORED ARRAY IS TRANSFORMED TO A VECTOR STORED SRiF ARRAY - 


RIN( (N+l)*(N+2)/2) 
ROUT( {N+i)*(N+2)/2) 


N 


^IS = 


IS = 


N.GT.O 

N.LT.O 


input Vector stored srif or u-n array 

OUTPUT IS THE CORRESPONDING U-C OR SRiF 
array (RIN=R0UT Is permitted) 

ABS(N)= matrix dimension 

THE (INPUT) SRIF ARRAY IS OUTPUT IN- U-D FORM 
THE (INPUT) U-D ARRAY IS OUTPUT IN SRiF FORM 
there Is NO RT. SIDE OR ESTIMATF. STORED IN 
column N+1 # AND RIN NEED HAVE ONLY. 

N columns. I.E. RlN(N*(N+l)/2) 

THERE IS A RT. SIDE INPUT TO THE SRIF AND 
AN ESTIMATE FOR THE U-D ARRAY. THESE RESIDE 
IN COLUMN N-fl. 


THIS SUBROUTINE USES SUBROUTJnF UTINV 

COGIZANT PERSONS G. J.BIERMAN/M.W.NEAD (JPL. FEB. 1977) 


RMDROOlO 

RUDR0020 

RUDR0030 

RUDROOAO 

RUDROOSn 

R1JDR0060 

RUOR0070 

RUDR0080 

RUOR0090 

RUDROlOO 

RUDROllO 

RUDR0120 

PUDR0130 

RUDR0140 

RUDR0150 

RUDR0160 

RUDR017O 

RUDR0180 

RUDROlOO 

RUDR0200 

RUDR0210 


IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION RIN(D# ROUT(l) 

ONE= 1.0 

NP1= IS + ABS(N) 

JJ=1 (3 INITIALIZE DIAGONAL INDEX 

IDIMR= NP1+(NP1 +l)/2 
IF (IS.EQ.l) RIN(IDIMR)= “ ONF 

IF (N.LT.O) GO To 30 
CALL UTIMV(RIN#NP1#R0UT) 

ROUT(l)= R0UT(1)**2 
IF (N.EO.l) return 
DO 20 J=2#N 


RUDR0220 

RUDR0230 

RUDR0240 

R1JDR0250 

PUDR0260 

RUDR0270 

PIJDR0280 

RUDR0290 

PUDR0300 

RUDR0310 

RDDP0320 

RUDR0330 

RUDR034Q 

RUDR0350 

RUDR0360 


S=ONE/ROUT ( JJ+J) 
ROUT(JJ+J)= ROUT ( JJ+J) **2 
JM1=J-1 
DO 10 I=1#JM1 

10 ROUT(JJ+I)= ROUT(JJ+n*S 
20 JJ=JJ+ J 
RETURN 


30 N=-N 

ROUT(l)= SQRT(RlN(l) ) 
IF(N.EO.I) GO TO 60 
DO 50 J=2»N 

ROUT(JJ+J)= SQRT(RINCJJ+J) ) 
S=ROUT (JJ+J) 

JM1=J-1 
DO 40 I=1#JM1 

40 ROUT(JJ+I)= RIN(JJ+I)*S 
50 JJ=JJ+J 

60 call UTINV(R0UT#NP1# ROUT) 


IDOGIBO-ITY OF THB 

Xpage b poor 


RETURN 

end 


RUDR0370 

RUDR03R0 

RUDR039n 

RUDR04n0 

PUDR0410 

RIJDR0420 

RUDR0430 

RUDR0440 

RUDR0450 

RUDR0460 

RUDR0470 

RUDR0480 

PUDR0490 

RUDR0500 

RUDR0510 

RUDR0520 

RUDR0530 

RUDR0540 

RUDR0550 

RUDR0560 

RUDR0570 

RUDR0580 

RUDR0590 
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SUBROUTINE THH(R»N»A»1A»M» SOS . MSTRT ) 

THHOOOlO 

THIS SUBROUTINE PERFORMS A DOUBLE PRECISION TRIANGULARIZATION THH00020 
OF A rectangular MATRIX INTO A SINGLY-SUBSCRIPTED ARRAY BY ‘ THH00030 

application of householder orthonormal transformations. THH00040 

THH00050 

R(N*(N+3)/2) VECTOR STORED SQUARE ROOT INFORMATION MATRIX THH00060 

(LAST N LOCATIONS MAY CONTAIN A RIGHT HAND SIDE) THH00070 
N NUMBER OF PARAMETERS THH00080 

A(IA»N+1) MEASUREMENT MATRIX THH00090 

lA Row DIMENSION OF A ThHOOIOO 

M NUMBER OF OBSERVATIONS IN THIS BATCH THHOOllO 

SOS accumulated Sum of squares of the RESIDUALS THH00120 

(Z-A*X(EST)**2) , INCLUDES A PRIORI THHD0130 

NSTRT FIRST COL OF THE INPUT A MATRIX THAT HAS A NONZERO THH00140 

entry, if nstrt.le.1i it is set to 1, This option thhooiso 

IS CONVENIENT WHEN PACKING A PRIORI RY BATCHES AND ThHOOIGO 
The a MATRIX HAS LEADING COLUMNS OF ZEROS. THH00170 

THHOOIBO 

ON ENTRY R CONTAINS A PRIORI SQUARE ROOT INFORMATION FILTER (SRIF) THHOOl 90 
ARRAY# AND ON EXIT IT CONTAINS THE A POSTERIORI (PACKED) ARRAY. THH00200 
ON entry a CONTAINS OBSERVATIONS WHICH ARE DESTROYED BY THE THH00210 

INTERhfAL COMPUTATIONS. THH00220 

ON entry if SOS IS .LT. ZERO .PROGRAM WILL ASSUME THERE IS NO THH00230 

RIGHT HAND SIdE DATA AND WILL NOT COMPUTF 505 OR USE LAST N THH00240 
LOCATIONS OF VECTOR R. THH00250 

THH00260 

COGNIZANT PERSONS G . J.RIFRMAN/N.hAMATa (JPLr OCT.1975) THH00270 

THH002RO 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) THH00290 

DIMENSION A(IA#1)#R(1) THH00300 

DOUBLE PRECISION SUM THH00310 

data ZERO/O.DO/» ONE/l.DO/ THH00320 

THHn0330 

IF ( NSTRT. LE.O) NSTRT=1 THH00340 

NP1=N+1 .Q no. columns of R THHD03S0 

IF(S0S,LT.ZER0) NP1=N Q NO COLS. = N IF SOS.LT.O THH00360 

KK=NSTRT*(N5TRT~1)/2 THH00370 

DO 100 J=NSTRT»N Q J-TH STEP OF HOUSEHOLDFR REDUCTION THHn0380 

KK=KK+J THH00390 

SUM=ZERO THH00400 

DO 20- 1 = 1 fM THHP0410 

SUM=SUM+A(I # J)+*2 THHn0420 

IF(SIJM.LE.ZERO) GO TO 100 Q IF J-TH CoL. OF A.EQ.O GO TO STEP J+lTMHn0430 
SUM=SUM+R (KK)**2 THH00440 

SUM=DSQRT(SUM) THH00450 

IF(R(KK) .GT.7ER0) SUMr-SUM THH00460 

DELTA=R(KK)-5UM THH00470 

R(KK)=SUH THH004B0 

BETA=ON£/(SUM*nELTA) THH00490 

JU=KK THH00500 

L=U THH00510 

J1=J+1 THH00520 


ON entry a CONTAINS OBSERVATIONS WHICH ARE DESTROYED BY THE 
INTERhfAL COMPUTATIONS. 

ON entry if SOS IS .LT. ZERO .PROGRAM WILL ASSUME THERE IS NO 
RIGHT HAND SIdE DATA AND WILL NOT COMPUTF 505 OR USE LAST N 
LOCATIONS OF VECTOR R. 

COGNIZANT persons G . J.RIFRMAN/N.hAMATa (JPLr OCT.1975) 

IMPLICIT DOUBLE PRECISION (A-H.O-Z) 

DIMENSION AdA.D.Rd) 

DOUBLE PRECISION SUM 
data ZERO/O.DO/r ONE/1. DO/ 


IF (NSTRT. LE.O) NSTRT=1 
NP1=N+1 

lF(SOS,LT.ZERO) NP1=N 
KK=NSTRT* ( NSTRT-l ) /2 
DO 100 J=NSTRT.N 
KK=KK+J 
SUM=ZERO 
DO 20- I = lrM 
20 SUM=SUM+A(I . J)+*2 

IF(SIJM.LE.ZERO) GO TO 100 
SUM=SUM+R (KK)**2 
SUM=DSQRT(SUM) 

IF(R(KK) .GT.7ER0) SUMr-SUM 

DELTA=R(KK)-5UM 

R(KK)=SUH 

BLTA=ON£/(SUM*nELTA) 

JJ=KK 

L=J 

J1=J+1 


13 NO. columns OF R 

ra NO COLS. = N IF SOS.LT.O 

Q J-TH STEP OF HOUSEHOLDFR REDUCTION 


+* READY TO APPLY J-TH HOUSEHOLDER TRANS. 
DO 40 K=JlrNPl 


THH00530 

THH00540 
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JJ=JJ+L 

THH00550 


L=L+1 

THH00560 


SUM=oELTA*R(JJ) 

THH00570 


DO 30 I=lrM 

THH00580 

30 

SUM=SUM+A ( I f J ) *A ( I r K ) 

THH00590 


IF(SUM.EQ.ZERO) go to 40 

THH00600 


SUM=SUM*BETA 

THH00610 


R ( J J ) =R { J J ) +SUM+DELT A 

THH00620 


DO 35 l=lfM 

THH00630 

35 

A(I»K)=A(I»K)+SUM*A(I.J) 

THH00640 

40 

CONTINUE 

THH00650 

100 

CONTINUE 

THH00660 


IF(SOS.LT.ZERO) RETURN 

THH00670 

THH006BO 


CALCULATE SOS 

THH00690 

TMH00700 


SUM=ZERO 

THH00710 


DO 110 I=1»M 

THH00720 

110 

SUM=SUM+A ( I » NPl ) **2 

THH00730 


SOS=DSQRT ( S0S++2+SUM ) 

THH00740 

THHn0750 


RETURN 

THH00760 


END 

THH00770 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


subroutine TRIMAT ( A,N.CAR|TExT,NCHAR,NAMES) 

TO display a Vector stored Upper triangular matrix in a 

TWO-DIMENSIONAL TRIANGULAR FqRMAT 

A(N*(n+i>/2) Vector containing upper triangular hajrix (Opj 

N dimension of matrix ( I j 

CaRIN) parameter names (I) 

TeXT( ) AN array of FielOaTA CHARACTERS TO BE PRINTED AS 

A title preceding THE MATRIX 

NcHAR number OF characters, including SPACES, that 

are to be printed in TEXT! ) 

ABs ( NCH aR) «LE, l26, NCHAR NEGATIVE IS USeD 

to avoid skipping to a new page to start 
printing 

names true to print Parameter names 

cognizant persons; 6. J.B iErMaN/M.W.NEaD (JPL, 0Ct»1V75) 

Double precision a(N) 

integer CaR(N), TEXT(I), U7), L1ST(7> 

logical names 

integer V(R),VFMT(7) 

data V/' (2X,* A6. IX.' '.'Dl7.6)V, 

1 VFMT/'7','Ol7X,6'.'o3‘»X,5','o5lX,R','06BX,3*,»oa5X,2»,M02X,I'/ 


Ml ,M2 

ROW 

LIMITS FOR EAcH 

PRINT sequence 

Ni ,M2 

COL 

limits FOR EAcH 

line of PRINT 

L{ I) 

LOC 

OF Each coLUM^g 

IN A ROW 

KT 

ROW 

counter 


KP 

PRINT counter 



• • * ♦ • initialize counters 

Ml = l 
M2 = 7 
Nlej 
KT = 0 
KPcq 

IF (.NOT, NAMES) V(2)hM5,2X' 

NC=IaBS(NCHAR)/6 
IF (M0D(NCHaR, 6) ,NE.D) NC'NC+i 
IF (nCHAR.Ge.O) write (6,200) ( TEXT ( I ) , I = 1 , NC ) 
IF (nCHaR^Lt^O) write (6,205) ( TEXT ( I ) , I = I , NC ) 

ID IF (M2.GT.N) M2=N 

IF (.NOT, NAMES) GO TO 20 
write (6,2lQ) (CAR( I ) , ItaNl ,M2) 

GO TO ‘(0 
20 M=N1 

L2=H2-N1+1 
DO 30 lBl,L2 
L1ST( I )«iM 
30 M=M+1 

WRITE (6,220) ( L I ST ( I) , 1 = 1 i L2 ) 


TRIMOOI D 
TRIH0Q2O 
TRIM0030 
TRIMOOHO 
TRIM0050 
TRIM0060 
TR'IMD070 

trimoobo 

TR1H0090 
TRIMOIOO 
TRiMOl 10 
TR1M0120 
TRIM0130 
TRIHOIRO 
TRIMOISO 
TRIM0I60 
TR 1M0170 
TRlMDiao 
TRIM0190 
TRIM0200 
TR IM0210 
TRIH0220 
TRIMD230 
TRIMD290 
TRIMD2S0 
TRIM0260 
TR IM0270 
TRIM0280 
TRIHQ290 
TR IM0300 
TRIH0310 
TRIM0320 
TRIM0330 
TR IH03H0 
TR IH035D 
TRIM0360 
TRIM0370 
TRIM03B0 
TRIM0390 
TRiMO'lOO 
TRiMUNlO 
TRIM0‘(20 
TRIMD‘»30 
TRIM04‘(D 
TRIM0'(5 d 
TRIH0‘(6d 
TR IHDM7o 
TR1M0‘(8d 
TR lM0*(9O 
TRiMOBOO 
TRIM051 0 
TRIM0520 
TRIM0530 
TRIM05'(0 
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‘♦0 CONTINUE 

••***• 

DO 190 IC«HI,H2 
K«l 

IF ( iciuEl (KT*7) » 60 TO 60 
iJiJoO 

DO 50 J«l,lC 
50 JJBJJ+J 

L{K)«UJ 
I 1«IC-Kt*7 

IF m.eaiTJ so jo 90 
60 TO 70 

60 continue ' 

I 1b1 

L(KJ-UK)f I 
70 continue 

DO 80 loll, 6 
K-K+1 
II«I+KT*7 

80 l(Kj«l{k-i )*n & obtain col index For row 

70 continue 

I2«MINO(8, (M2+1-KT*7) )-I I 
V(3)“VFMT III) 

IF I.NOt.NAHES) 60 TO 180 

WRITE (6»V) CARnC) , UtU I ) ) ,i»i , 12) 

GO TO I90 

lao WRITE 16,V) IC,(A(L( I >) I I»I , l2) 

170 continue 

IF (M2.EQ,N) return 
N 1 0 M 2 + 1 
M2*>M? + 7 
KT*KT+1 

Kp*>Kp + 1 

IF (KP.LT.3) 60 TO Iq 

write (6,200) (TEXt( I ) , I« l ,NC j 

60 TO 10 


200 

format 

(1H1,2X,21A6) 

0 TITLE 

205 

format 

{1H0,2X,2S A6) 

0 TITLE 

210 

format 

(IH0,5X,7UiX,A6)) 

0 horizontal NaMeS 

220 

format 

end 

(1H0,3X,7( 1 IX, 16) ) 



TRIM05S0 

TR1M0560 

TRJM057Q 

TR1m0580 

TRIH0S9O 

TRIM0600 

TRlM06i0 

TR IM0620 

TRIH0630 

TRlH06‘(0 

TRIH0650 

TRIM0660 

TR IH0670 

TRIHO 680 

TRIM0670 

TRIM0700 

TRIM0710 

TRIM0720 

TRIH0730 

TR lM07*t0 

TRIM0750 

TRIM0760 

TRIM0770 

TR1H0780 

TRIM0790 

TRIH0800 

TRIM0810 

TRIHO820 

TRIM083Q 

TR iMOe'iO 

TRIH0850 

TRIH0860 

TRIM0870 

TRlMOBSO 

TRIH0B90 

TRlM090D 

TRIMD910 

TRIM0920 

TRIM0930 

TRIH09‘)0 

TR IMD95D 

TR IM0960 

TRIM097D 

TRIH0980 

TRIH0970 
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SUBROUTINE TTHH(RrRA»N) 

■ This SUBROUTINE COMBINES TWO SINGLE SUBSCRIPTED SRIF ARRAYS 
USING HOUSEHOLDER ORTHOGONAL TRANSFORMATIONS 

R{N*(N+l)/2) VECTOR STORED SRIF ARRAY. 

RA(N+(N+l)/2) THE SECOND VECTOR STORED ‘5RIF ARRAY 
N DIMENSION OF THE ESTIMATED PARAMETER VECTOR. 

A NEGATIVE VALUE FOR N IS USED TO NOTE THAT 
R AND RA Have RT. hand sides included AND 
HAVE DIM=ArsN*CABSN+3)/2. 

ON EXIT RA IS CHANGED AND R CONTAINS THE RESULTING SRIF ARRAY 

COGNIZANT Persons G. J.BIERMAN/M.W.NEAD UPL. JaN.1976) 

IMPLICIT double PRECISlON(A-HrO-Z) 
dimension RA(1) . R(1 ) 

double precision sum q For use in single precision version 

ZERO=0 . 

ONE=l. 

NP1=N 

IF (N.GT.O) go To 10 

N=-N 

NP1=N+1 

10 1JS=1 D IJCSTART) 

KK=0 

DO 100 J=1»N (S J-Th step of HOUSEHOLDER REDUCTION 

KK=KK+J 
SUM=R(KK)**2 
DO 20 I=IJS»KK 
20 SUM=SUM+RA ( I ) **2 

IF (SUM. LE. ZERO) GO TO lOO 
SUM=SQRT(SUM) 

IF (R(KK) .GT.ZERO) SUM=-SUM 

DELTA=R(KK)-5UM 

R(KK)=SUM 

bETA=ONE/ (SUM*DELTA) 

JJ=KK 

L=J 

UP1=J+1 

IKS=KK+1 

* * * J-TH HOUSEHOLDER TRANS. DEFINED 

40 LOOP APPLIES TRANSFORM. TO COLS. J+1 TO NPl 
DO 40 K=JPl.NPi 
JJ=JJ+L 
L=L+1 
IK=IKS 

sum=delta*r(jj) 

DO 30 I=lJSfKK 
SUM=SUM+KA ( IK ) *RA ( I ) 

30 IK=IK+1 

IF (SUM. EQ. ZERO) GO TO 40 
SUM=SUM*BETA 


TTHHOOlO 

TTHH0020 

TTHH0030 

TTHH0040 

TTHH0050 

TTHH0060 

TTHH0070 

TTHH0080 

TTHH0090 

tthhoioo 

TTHHOIIO 

TTHH0120 

TTHH0130 

TTHH0140 

TTHH0150 

TTHH0160 

TTHH0170 

TTHHOIBO 

TTHH0190 

TTHH0200 

TTHH0210 

TTHH0220 

TTHH0230 

TTHH0240 

TTHH0250 

TTHH0260 

TTHH0270 

TTHH0280 

TTHH0290 

TTHH0300 

TTHH0310 

TTHH032Q 

TTHH0330 

TTHH0340 

TTHH0350 

TTHH0360 

TTHH0370 

TTHH0380 

TTHH0390 

TTHH0400 

TTHH0410 

TTHH0420 

TTHH0430 

TTHH0440 

TTHH0450 

TTHH0460 

TTHH0470 

TTMH0480 

TTHH0490 

TTHH0500 

TTHM0510 

TTHH0520 

TTHH0530 

TTHH0540 


83 



77-26 



R { J J ) =R { JJ ) +SUM*DELTA 

TTHH0550 


IK=IKS 

TTHH0560 


DO 35 I=IJS»KK 

TTHH0570 


RA ( IK ) =RA ( IK) +SUM+RA ( I ) 

TTHH05B0 

35 

IK=IK+1 

TTHH0590 

40 

IKS=IKS+K 

TTHH0600 

100 

IJS=KK+1 

TTHH0610 

TTHH0620 


RETURN 

TTHH0630 


END 

TTHH0640 
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SUBROUTINE T?ERO (RrN»IS,lF) 

T^EROOIO 

TO ZERO OUT Rows IS CiSTART) TO IF ClFiNAL) OF A VECTOR TZER0020 

STORED UPPER TRIANGULAR MATRIX T7ER0030 

T7ER0040 

R(N*(N+1)/2) input Vector stored upper triangular matrix tzerooso 

N DIMENSION Op R T7ER0060 

IS FIRST ROW OF R THAT IS TO RE SET TO ZERO T7ER0070 

IR LAST ROW OF R THAT IS To BF SET TO ZERO T7ER0080 

T7ER0090 

COGNIZANT persons: G.J.BIERMAN/C.F, PETERS (JPL» NOV. 1975) T7ER0100 

T7ER0110 

IMPLICIT DOUBLE PRECISION (A-H,0-Z) TZER0120 


DIMENSION R(l) T7ER0130 

C TZEP0140 

2£R0=0.Q T7ER0150 

IJS=IS*(IS-1)/? T7ER0160 

DO 10 I=IS»IF T7ER0170 

IJS-US+I T7ER0180 

IJ=IJS T7ER0190 

DO 10 J=IrN TZER0200 

R(IJ)=ZER0 T7ER0210 

1J;=IJ+J TZER0220 

10 CONTINUE T7ER0230 

C T7ER0240 

RETURN T7FR0250 

end T7ER0260 
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C 


C 


SUBROUTINE UDMES {U»NrRfA»G» ALPHA) 

COMPUTES estimate AND U-D MEASUREMENT UPDATED 
COVARIANCE* P=UDU**T 

*♦* INPUTS *** 

U UPPER TRIANGULAR MATRIX# WITH D ELEMENTS STORED AS THE 

diagonal. U IS vector STORED AND CORRESPONDS TO THE 
A PRIORI COVARIANCE. IF STATE ESTIMATES ARE COMPUTED# 
THE LAST COLUMN OF U CONTAINS X, 

N DIMENSION OF THE STATE ESTIMATE. 

R MEASUREMENT VARIANCE 

A VECTOR OF MEASUREMENT COEFFICIENTS# IF DATA THEN A(N+1 

ALPHA IF ALPHA LESS THAN ZERO NO ESTIMATES ARE COMPUTED 
(AND X AND Z NEED NOT BE INCLUDED) 

OUTPUTS *** 

U UPDATED# VECTOR STORED FACTORS AND ESTIMATE AND 

U( (N+1) (N+2)/2) contains (Z-A**T*X) 

ALPHA INNOVATIONS VARIANCE OF THE MEASUREMENT RESIDUAL 
G VECTOR OF UNWEIGHTED KALMAN GAINS# K=G/ALPHA 

A CONTAINS U**TA AND ( Z-A**T*X ) /ALPHA 

COGNIZANT Persons: g.j. bierman/m.w. nead (jpl# SEPT.iqye) 

IMPLICIT DOUBLE PRECISION (A-H#0-Z) 

DIMENSION U(l)f A(i)# G(l) 

DOUBLE PRECISION SUM 
LOGICAL lEST 

ZERO=0.0 


udmesoio 

UDMES020 

UDMES'030 

UDMESQ40 

UDMES050 

UDMES060 

UDMES070 

UDMES080 

UDMES090 

UDMESiOO 

UDMESllO 

UDMES120 

UDME5130 

)=ZUDMES140 

UDMES150 

UDMES160 

UDMES170 

UDMES180 

UDMES190 

UDMES200 

UDMES210 

UDMES220 

UDME5230 

UDMES240 

UOMES250 

UDMES260 

UDMES270 

UOMES280 

UDMES290 

UDMES300 

UDMES310 

UDMES320 

UDMES330 

UDMES340 


l£5T=. FALSE. 

OUE=i. 

NP1=N+1 

NT0T=N*NP1/2 

IF (ALPHA.lt. ZERO) GO TO 3 
SUM=A(NP1) 

DO 1 J=i#N 

1 SUM=SUM-A(J)*U(NT0T+J) 

U(NTOT+NPl)=SUM 0 Z=Z-A+*T*X 

l£ST=.TRUE, 

3 KJ=NT0T 

DO. 10 J=N»2#-1 
SUM=A(J) 

UM1=J-1 

DO 5 K=JMl,i#“l 
KJ=KU-1 

5 SUM=SUM+U(KJ)*A(K) 

A(J)=SUM 

KJ=KJ-1 

10 G(J)=SUM*U(KJ+J) 


UOMES350 

UDMES360 

UDMES370 

UDMES380 

UDMES390 

UDMES400 

UDMES410 

UDMES420 

UDMES430 

UDMES440 

UOMES450 

UOMES460 

UDMES470 

UOMES480 

UDMES490 

UDMES500 

UDMES510 

UDMES520 

UDMES530 

UDMES540 

UDMES550 
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G(l)=Ua)*A(l) 

A=U**T*A and 6=D*{U**T*A) 

SUM=R+G ( 1 ) *=A ( 1 ) 

gamma=o 

IF (GCl) .EG,.ZERO) GO TO U 

gamma=one/sum 

Ud)=Ua)*R*GAMMA 

11 KJ=2 

DO 20 J=2»N 
BETA=SUM 

SUM=SUf/|+G{ J)*A(J) 
P=-A(J)+GAMMA 
JM1=J-1 
DO 15 K=1»JM1 
S=U(KJ) 

U(KJ)=S+P»G(K) 

G(K)=G (K)+G( J)*S 
15 KJ=KJ+1 

IF (G(J) .EQ.ZERO) GO TO 20 
GAMHA=0NE/SUm 
U ( K J ) =U ( K J ) +BETA*6AMMA 
20 KJ=KJ+1 
ALPHA=SUM 


D SUm(1) 

0 FOR R=0 CASE 
Q FOR R=0 CASE 

D D(l) 


ra BETA=SUiV(J-1) 

9 SUm(J) 

9 P=-F(J)*(1/SUM( J-1) ) EQn(2D 


0 EON (22) 

0 EON (23) 

0 FOR R=0 CASE 
0 GAMMA=1/SUM(J) 

0 D(J) EQN(19) 


paper r PP. 337-346, 


EQN. NOS. REFER TO BIERMAN’S 1975 CDC 


IF (.NOT.iEST) Return 
A (NP l ) =U (NTOT+ nPI ) +GAMMA 
DO 30 J=lfN 

30 U (NT0T+J)=U (nTOT+J ) +G ( j ) *A (NP l ) 
C 

RETURN 

END 


UDMES560 

UDMES570 

UDMES580 

U0MES590 

UDME5600 

UDI^ESeiO 

UOMES620 

UDMES630 

UDMES640 

UDMES650 

UDMES660 

UDMES670 

UDME5680 

UDMES690 

UDMES700 

UDMES710 

UDMES720 

UDMES730 

UDME5740 

UDMES750 

UDMES760 

UOMES770 

UOMES780 

UDMES790 

UOMES800 

UOMES810 

UDMES820 

UDMES830 

UDMES840 

UDMES850 

UDMES860 

UDHES870 

UDMES880 

UOMES890 

UDMES900 


87 



onoooooooooon 


77-26 


C 


c 

c 


c 


subroutine UD2C0V (UINrPOUT»N) 

TO obtain a covariance from its U-D FACT0RI7ATI0N. ROTH MATRICES 
ARE VECTOR STORED AND THE OUTPUT COVARIANCE CAN OVERWRITE THE 
INPUT U-D ARRAY. UIN=U-D IS RELATED TO POUT VIA POUT=UDU{**T) 

UIN(N*(N+1)/?) INPUT U-D FACTORS. VECTOR STORED WITH THE D 
entries stored ON THE DIAGONAL OF UIN 
P0UT(N*(N+l)/2) OUTPUT COVARIANCE. VECTOR STORED. 

(POUT=UIN IS PERMITTED) 

N DIMENSION OF THE MATRICES INVOLVED 

COGNIZANT persons: 6. J ,BIERMan/M , W. NEAD (JPL» FEB. 1977) 

IMPLICIT DOUBLE PRECISION tA-H.O-Z) 

DIMENSION UIN(l). POUT ( 1 ) 

POUT(l)=UIN(l) 

JJ=1 

DO 20 J=2.N 

JJL— <JU Q (J— l.J— l) 

U J— J 

P0UT(JJ)=UIN(JJ) 

S=POUT( JJ) 

11=0 

JM1=J-1 
DO 20 1=1. JMl 
II=II+I 

ALPHA=S*UIN(JJL+I) 13 JJL+I=(I.J) 

IK=II 

DO 10 K=I.JM1 

POUT(IK)=POUT(IK)+aLPHa*UIN(JJL+K) (3 JJL+K=(K,J) 

10 IK=IK+K 

20 POUT(JJL+I)=ALpHA 

RETURN 

end 


UD2C0010 

UD2C0020 

UD2C0030 

UD2C0040 

UD2C0050 

UD2C0060 

UD2C0070 

UD2C0080 

UD2C0090 

Un2COiOO 

UD2C0110 

UD2C0120 

UD2C0130 

UD2CO140 

UD2C0150 

UD2C0160 

UD2C0170 

UD2C0180 

UD2C0190 

UD2C0200 

UD2C021D 

UD2CO220 

UD2C0230 

UD2C0240 

U02C0250 

UD2C0260 

UD2COP70 

UD2C0280 

UD2C0290 

UD2C030Q 

UD2C0310 

UD2C0320 

UD2C0330 

Un2C0340 

UD2CO3S0 

UD2C03ftO 

UD2C0370 

UD2C0380 
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subroutine UD2SlG(UfN»SlGrTEXT,NCT) 

compute standard deviations (STSMAS) from u-d covariance factors 

U(N*(N-H)/2) INPUT vector STORED ARRAY CONTAINING THE U-D 

factors. The d (diagonal) elements are stored 
ON THE diagonal 

SIG(N) VECTOR OF OUTPUT STANDARD DEVIATIONS 

TEXT( ) ARRAY OF FIELDATA CHARACTERS TO BE PRINTED 

PRECEDING THE VECTOR OF SIGMAS 
NCT number of CHARACTERS IN TEXTf 0.LE.NCT.LE.126 

IF NCT=oi NO Sigmas are printed 

COGNIZANT persons: G. J.BIERMaN/M. W, NEAD {JPL. FEB. 1977) 

IMPLICIT DOUBLE PRECISION (A-H*0-2) 

INTEGER TEXT(l) 

DIMENSION U(l)r SIG(I) 

JJ=1 

SIG(1)=U(1) 

DO 10 J=2»N 

UUL— UU Q (U"*1^U— 1) 

S=U(JJ) 

SIG(J)=S 

JM1=J-1 

DO 10 I=1»JM1 

10 SIG(I)=SIG(I)+S*IJ{JJL+I)**2 
WE NOW HAVE VARIANCES 
DO 20 J=1»N 

20 SIG(J)=SQRT(SIG(J) ) 

IF (NCT.EQ.O) go TO 30 
NC=NCT/6 

IF (MOD(nC»6) .NE.O) NC=NC+1 
write (6»40) (TEXTCI) »I=1»NC) 

WRITE (6f50) (SIGH) f I=1»N) 

30 RETURN 

40 format (lH0r2X.2lA6) 

50 FORMAT (IHOr (6D18.10) ) 

END 


UD2SI010 

UD2SI020 

UD2SI030 

UD2SI040 

UD2SI050 

UD2SI060 

Un2SI070 

UD2SI080 

UD2S1090 

UD2SI100 

UD2SI110 

UD2SI120 

UD2SI130 

UD2SI140 

UD2SI150 

UD2SI160 

UD2SI170 

UD2SI180 

UD2SI190 

UD2SI200 

UD25I210 

UD2SI220 

UD2SI230 

UD2SI240 

UD2SI250 

UD2SI260 

Un25I270 

Un2SI280 

UD2SI290 

UD2SI300 

UD2SI310 

UD2SI320 

UD2SI330 

Ud2SI340 

UD2SI350 

UD2S1360 

UD2SI37Q 

UD2SI380 

UD2SI390 

UD2SI400 

UD2SI410 

Un2SI420 

UD2SI430 

UD2SI440 
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SUBROUTINE UTINV (RIN»NrROUT) 

TO INVERT AN UPPER TRIANGULAR VECTOR STORED MATRIX AND STORE 
the RESULT IN VECTOR FORM. THE ALGORITHM IS SO ARRANGED THAT 
THE RESULT CAN OVERWRITE THE INPUT. 

IN addition To solve RX=Z» set RINtN*(N+l)/2+l)=Z(l) * ETC.» 
AND SET RIN{ (N+1)*(N+2)/2)=-1. CALL THE SUBROUTINE USING N+1 
INSTEAD OF N. ON RETURN THE FIRST N ENTRIES OF COLUMN N+1 
WILL CONTAIN X. 

RIN(N*(N+l)/2) INPUT VECTol? STORED UPPER TRIANGULAR MATRIX 
N MATRIX DIMENSION 

ROUT (N* (N+1 )/2) OUTPUT VECTOR STORED UPPER TRIANGULAR MATRIX 
inverse 

COGNIZANT PERSONS! G. J.BTERMAN/J.ELLIS (JPL» SEPT. 1976) 

double PRECISION RIM(l)r ROtJT(l)f WORK» ONF» ZERO .DIN 
DATA ONE/1. ODO/.ZERO/ O.ODO/ 

IPV = N»(N+l)/2 

IN = IPV 
DO 6 1=1. N 

IF (RIN(IPV) .NE.ZERO) go TO 1 

WRITE (6.10) I 

RETURN 

1 DIN = One/ RIN(IPV) 

ROUTC IPV ) = DIN 

MIN =N 
KEND = 1-1 
LANF = N - KEND 

IF (I.EQ.l) GO TO 5 

2 J= IN 

INITIALIZE ROW LOOP 

DO 4 K=1.KEND 
WORK =ZERO 
MIN= MIN - 1 

LIN= IPV 
LOT= J 

START INNER LOOP 

DO 3 L=LANF. MIN 
LI,N= LIN+L 
LOT= LOT+1 

3 WORK = WORK + RlN(LIN)* ROUT(LOT) 

ROUT(U) = -WORK* DIN 

4 J= J- MIN 

5 IPV = IPV -MIN 

6 IN= IN -1 
RETURN 

10 format (IHO. lOX. 'UTINV diagonal’ . 14. » IS ZERO’) 

END 


UTINVOlO 
UTIMV02Q 
UTINV030 
UTINV040 
UTINV050 
UTINV060 
UTINV070 
UTINVOSO 
UTINV090 
UTINVIOD 
UTINVllO 
UTINV120 
UTTNV130 
UTINV140 
UTINV150 
UTINV160 
UTINV170 
UTINVISO 
UTINV190 
UTINV200 
UT1MV210 
UTINV220 
UTINV230 
UTINV240 
UTINV250 
UTINV260 
UTINV270 
UTINV280 
UTINV290 
UTINV300 
UTINV310 
UTINV320 
UTINV330 
UTINV340 
UTINV350 
UTIMV360 
UTINV370 
UT1NV380 
UTINV 390 
UTINV400 
UTINV410 
UTINV420 
UTINV430 
UT1NV440 
UTINV450 
UTIMV460 
UTINV470 
UTINV4B0 
UTINV490 
UTINV500 
UTINV510 
UTINV520 
UTINV530 
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C 

.SUBROUTINE IITIROW (RINrNrROUT»NRY) 


TO COMPUTE THE INVERSE OF AN UPPER TRIANGULAR (VECTOR STOREn) 
MATRIX WHEN THE LOWER PORTION OF THE INVERSE IS GIVFN 

ON INPUT: ' 


RX 

RiNr 

RXY 

♦ 

R0UT= 

♦ 

where 

RX 

R= 

RXY 

* 

♦ 

0 

RY**-1 

^ 0 

RY 


ON OUTPUT: RIN IS UNCHANGER AND ROuT=R+*-l 
» THE RESULT CAN OVER-WRITF THE INPUT (I.E. RIN=R0UT) 


RIN(N*(N+1)/?) 

N 

R0UT(N*(N+l)/2) 


NRY 


INPUT vector STORED TRIANGULAR MATRIX 
THE BOTTOM NRY ROWS ARE IGNORED 
MATRIX dimension 

OUTPUT VECTOR STORED MATRIX. ON INPUT THE 
BOTTOM NRY ROWS CONTAIN THE LOWER PORTION 
OF R**~l. ON OUTPUT R0UT=R*+-1 
DIMENSION OF LOWER (ALREADY INVERTED) 
TRIANGULAR R. IF NRY=0, ORDINARY MATRIX 
INVERSION results. 


cognizant persons: G.J.BIERMAM/M.W.NEAD (JPL MARCH 1977) 


DOUBLE PRECISION RIN(1)» ROUT(l)» SUM* ZERO* ONE* DiNV 
DATA ONE/1. DO/* ZERO/O.DO/ 


INITIALIZATION 


10 


C 


C 



NK=N*(N+1)/2 
ISTRT=N“NRY 
IHLST=ISTRT+1 
Il=ISTRT*IRLST/2 
DO 40 IR0W=ISTRT, 1 *-1 

IF (RIN(II) .NE.ZERO) GO 
viiRlTE (6*50) IROW 
RETURN 

DINV=ONE/RIN(Il) 

ROUT(II)=DTNV 

KJS=NR+IROW 

IKS=II+IROW 


0 NO. elements in R 
0 FIRST ROW TO BE INVERTED 
N IRLST=PREVIOUS IROW INDEX 
U II=DIAGONAL 

J 0 


a KJ( START) 
Q IK (START) 


IF (IRLST.GT.N) GO TO 35 
DO 30 J=N*IRLST*-1 
KJS=KJS“J 
SUM=ZERO 
IK=IKS 
KU— K J5 


DO 20 K=IRLST*J 
KJ=KJ+1 

5UM=SUM+RIN(IK)*R0UT(KJ) 


UTinoolo 

UTIR0020 

UTIR0030 

UTIR0040 

UTIR0050 - 

UTIR0060 

UTIR0070 

UT1R0080 

UTIR0090 

UTIPOIOQ 

UTIROllO 

UTIR0120 

UTIR0130 

UTIR0140 

UTIR0150 

UTIR0160 

UTIR0170 

UTIROiaO 

UTIR0190 

UT1R0200 

UTIR0210 

UTIR0220 

UTIR0230 

UTIR0240 

UTIP0250 

UTIR0260 

UTIR0270 

UTIR02RO 

UTIR0290 

UTTR0300 

UT1R0310 

UTIR0320 

UTIR0330 

UTIR0340 

UTIR0350 

UTIR0360 

UTIRO370 

UTIPO380 

(ITIR0390 

UTIR0400 

UTIRO410 

UTIR0420 

UTIR0430 

UTIRO440 

UTIR0450 

UTIRO460 

UTIR0470 

UTIRO450 

UTIR0490 

UTIR0500 

UTIR0510 

UTIR0520 

UTIRO530 

UTIR0540 
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20 IK=IK+K UTIR05S0 

C UTIR0560 

30 ROUT(KJS)=-SuM*DINV UTIR0570 

35 IRLST=IROW UTIR0580 

40 II=II-IROW UTIR0590 

RETURN • UTIR0600 

50 format (IHOflOXURlN DIAGONAL* » I4r * IS ZERO*) UTIR0610 

END UTIR0620 
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C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


40 


50 


C 


C 

c 

c 


60 

70 

100 


110 

C 


SUBROUTINE WGS (W» iMAX'Wf IWVJw.DiUtV) . WGSOOOlO 

MODIFIED GRAMM-SCHMIDT ALGORITHM FOR REDUCING WDW(**T> TO UnU(«*T) WGS00020 
FORM WHERE U IS A VECTOR STORED TRIANGULAR MATRIX WITH THE WGS00030 

RESULTING D ELEMENTS STORE On THE DIAGONAL ' ' WGS00040 


w(iw*Jw) Input matrix to be reduced to triangular form. 

THIS MATRIX IS destroyed BY THE CALCULATION " 
IW.LE.IMAXW. " 

DtlW) VECTOR OF NON-NFGATIVE HEIGHTS FOR THE 

ORTHOGONALlZATlON PROCESS. THE D»S ARE UNCHANGED 

BY THF CALCULATION. 

U(IW*(IW+1)/?) OUTPUT UPPER TRIANGULAR VECTOR STORED OUTPUT 

V(JW) WORK vector 

(SEE rook: 

* factorization methods for discrete sequential estimation »» 

BY G.J.BIERman) 

estimation 

COGNIZANT persons: G. J.BIERMAN/M.W.NEAD (JPL» march 1977) 

IMPLICIT double precision (A-h*0-Z) 

DIMENSION W(IMAXWri)» D(l)» U(l)» V(l) 


Z=0.Q 


WGS00050 

WGS00060 

WGS00070 

WGS00080 

WGS00090 

WGSOOlOO 

WGSOOllO 

WGS00120 

WGSnoiao 

WGS00140 

WGsnoiso 

WGS00160 
WGS00170 
WGS005B0 
WGSn019Q 
WGS 00200 
WGSOOPIO 
WGS00220 
WGS00230 
WGS00240 
WGS00250 


ONE=l . 0 

00 100 J=1W»1»-1 
SUM=Z 

DO 40 K=lfJW 
V(K)=W(J»K) 
U(K)=D(K)*V(K) 
SUM=V(K)*UtK)+SUM 
W( J* J)=SUM 

IF (J.EQ.i) GO TO 100 
DINV=Z 

IF (SUM.GT.Z) nlNV=ONE/SUM 
JM1=J-1 


WGS00260 

WGS00270 

WGS00280 

WGS00290 

W6S00300 

QU HERE IS USED AS A WORK VECTOR WGS00310 

WGS00320 

[3,E0.(4.9) OF BOOK* NEW C(J) WGS00330 

WGS00340 

WGS00350 

WGS00360 

WGS00370 


DO 70 K=1»JM1 
SUM=2 


WGSn0380 

WGS00390 


DO 50 1=1 »JW 

SUM=W(K» D*U(I)+SUM 
SUM=SUM*DINV 


WGSn0400 

WGS00410 

W6S00420 


WG500430 


DO 60 I=1»JW 

W(K»I)=W(K,I)-SUM*V(I) 


WGS00440 

WGSnOUSO 


W(J*K)=SUM B EQ.(4.ln) OF BOOK 

CONTINUE a U(K»J) STORED IN W(J»K> 

THE LOWER PART OF W IS U TRANSPOSE 


U=o 

DO 110 J=1»IW 
DO 110 I=1*J 
IJ=IJ+1 
U(IJ)=W(J*I) 

RETURN 93 

END 


WGS00460 
WGSPQ470 
WGS 00 480 
WGSn0490 
WSS00500 
WGSOOSIO 
W6Sf)05P0 
WGSnosso 
WGS00S40 
WGSOO5S0 
WGSOOSGQ 
WGS00570 
WfiSn05fl0 
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